Multivariate statistics (classical methods)

Thomas Fuß/Gabriel Singer

WS24/25

Multivariate statistics - in R

What you have learned so far:

Ecological datasets have one or both of these datasets (if not more):

A short intro to classical ordination methods - in R

Ordination is an umbrella term for multivariate methods which arrange data along a gradient (axis, scale). A sloppy explanation is “putting things in order”.

Ordination is a word almost exclusively used by ecologists. Outside most of these methods are known as scaling techniques. The phrase gradient analysis makes a bit more sense in ecology.

Principal component analysis (PCA)

PCA can be understood as a MLR on theoretical (latent) instead of observed dependent variables.

The regression coefficients in PCA are factor loadings.

The predicted values of the theoretical dependent are scores.

For instance, in a PCA on this matrix of n objects times p variables PC1 is computed as:

\[ PC_1=f_1X_1+f_2X_2+...+f_pX_p \] There must be as many PCs as original variables. The factor loading matrix translates p \(X\)-variables into p \(PCs\).

PCA can be understood as a rotation of the original coordinate system made up of p \(X\)-variables into a new coordinate system defined by p \(PCs\).


Principal component analysis (PCA) in R

library(vegan)
## Lade nötiges Paket: permute
## Lade nötiges Paket: lattice
## This is vegan 2.6-4
library(shape) # plots nice arrows
setwd("Z:/COURSES/datenanalyse/vu_datenanalyse/unit_3") #set working directory
mara.raw<-read.table(file="data/MaraRiver.txt",header=TRUE) # water chemistry in 54 streams, 3 types of land use
dim(mara.raw)
## [1] 54 49

head(mara.raw)
##   site sampling.time diff.solar.noon landuse order     Q canopy Temp  cond   pH
## 1 Aa20         10:04          145.92       A     2 3.572    0.2 13.4 229.1 7.45
## 2 Aa21         14:01           91.08       A     1 6.600    0.2 17.0 205.4 8.03
## 3 Aa22         16:35          245.08       A     1 1.650    0.3 14.5 205.7 7.88
## 4 Aa23         10:05          145.22       A     1 5.472    0.4 12.1 181.8 7.96
## 5 Aa6b         17:03          270.78       A     2 1.800    0.2 17.5 404.4 7.79
## 6  Af5         14:04           94.38       F     3 4.200    0.4 14.5 129.4 7.94
##     turb  Alk       TSS      PO4         NH4      NO2         NO3       Br
## 1   0.02 2160  12.55102 5.16e-07 0.000022800 8.91e-07 0.000054500 2.45e-06
## 2  35.50 2160  59.44444 8.32e-07 0.000023600 6.08e-07 0.000007810 2.06e-06
## 3   0.02 2160  10.57143 1.37e-07 0.000027800 1.13e-06 0.000047300 1.39e-06
## 4   0.02 2160  26.96629 1.00e-06 0.000020700 2.17e-06 0.000124766 1.24e-06
## 5 388.00 3740 209.41176 1.37e-07 0.000446205 1.98e-05 0.000092700 2.47e-06
## 6  28.80 2160  10.21053 1.37e-07 0.000019300 6.95e-07 0.000034900 8.64e-07
##            Cl       Fl      SO4          Na           K          Ca          Mg
## 1 0.000265196 4.47e-05 8.58e-05 0.000930796 0.000423884 0.000328085 0.000152685
## 2 0.000237582 3.91e-05 5.09e-05 0.000843062 0.000323802 0.000323170 0.000146390
## 3 0.000220969 2.87e-05 3.56e-05 0.000622227 0.000373881 0.000372199 0.000157663
## 4 0.000215891 3.03e-05 4.66e-05 0.000716442 0.000299606 0.000300339 0.000132606
## 5 0.000354103 6.12e-05 5.97e-05 0.001070465 0.000989155 0.000507910 0.000253652
## 6 0.000124757 2.87e-05 5.13e-05          NA          NA          NA          NA
##        TDN  SUVA254      slope slope_short  SR_Helms E2.to.E3 E4.to.E6      FIX
## 1 7.82e-05 2.486708 0.01484335  0.01508883 0.8385724 5.729730 7.333333 1.432674
## 2 3.20e-05 2.749288 0.01304546  0.01313918 0.8853058 4.643519 5.000000 1.352617
## 3 7.62e-05 2.957983 0.01374024  0.01376633 0.8516180 4.986425 8.000000 1.371267
## 4 1.48e-04 3.069472 0.01360719  0.01317551 0.8079401 4.892857 6.714286 1.349814
## 5 5.59e-04 1.722983 0.01444961  0.01542405 0.8943684 5.288889 7.833333 1.481364
## 6 5.49e-05 3.088685 0.01376586  0.01382645 0.8333595 4.981043 6.428571 1.336503
##       FIX2       HIX beta.alpha      X7c_c1    X7c_c2    X7c_c3    X7c_c4
## 1 1.781212 0.8178391  0.7186987 0.002585346 0.5892004 0.4082030 0.2100508
## 2 1.709486 0.7345711  0.6978052 0.075614095 0.7089239 0.4878008 0.2858983
## 3 1.726269 0.7851819  0.6829999 0.037594156 0.8459288 0.5640788 0.3114274
## 4 1.701107 0.8066407  0.6676339 0.010937179 0.6222450 0.4279429 0.2482675
## 5 1.961382 0.6234517  0.7639094 0.222123376 1.8301573 1.7719178 0.5627895
## 6 1.722147 0.8536532  0.6777914 0.035057086 1.3009351 1.0052846 0.5806384
##       X7c_c5    X7c_c6    X7c_c7 humF.to.protF humF.per.DOC protF.per.DOC
## 1 0.03701974 0.3347390 0.1175206      9.815030    0.3153769    0.03213204
## 2 0.09037746 0.3795044 0.3062165      3.943447    0.2652603    0.06726611
## 3 0.05079516 0.4395259 0.2434638      6.511800    0.3026556    0.04647802
## 4 0.04149694 0.3223648 0.1405717      8.397779    0.2706781    0.03223210
## 5 0.63429053 1.1599349 1.5369801      2.224790    0.4052359    0.18214567
## 6 0.00000000 0.7400336 0.2134428     14.595144    0.5545706    0.03799693
##    DIC.d13C     pCO2    EpCO2
## 1 -2.673634       NA       NA
## 2 -8.590664 1184.457 3.060612
## 3 -4.565283 1753.501 4.531012
## 4 -7.214370 1094.980 2.829405
## 5 -4.717236 3730.961 9.640726
## 6 -9.946630       NA       NA
mara<-mara.raw[,c(4,9:26)] #only select columns of interest

which(is.na(mara), arr.ind=TRUE) #find rows with NAs
##      row col
## [1,]   6  15
## [2,]   6  16
## [3,]   6  17
## [4,]   6  18
mara<-mara[-6,] #delete row with NAs

landuse<-mara$landuse
wc<-mara[,-1]

apply(wc, 2, hist) #check normal distribution

## $cond
## $breaks
## [1]  50 100 150 200 250 300 350 400 450
## 
## $counts
## [1] 29 12  8  3  0  0  0  1
## 
## $density
## [1] 0.0109433962 0.0045283019 0.0030188679 0.0011320755 0.0000000000
## [6] 0.0000000000 0.0000000000 0.0003773585
## 
## $mids
## [1]  75 125 175 225 275 325 375 425
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $pH
## $breaks
## [1] 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8
## 
## $counts
## [1]  3 11 18 13  6  0  1  1
## 
## $density
## [1] 0.28301887 1.03773585 1.69811321 1.22641509 0.56603774 0.00000000 0.09433962
## [8] 0.09433962
## 
## $mids
## [1] 7.3 7.5 7.7 7.9 8.1 8.3 8.5 8.7
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $turb
## $breaks
## [1]    0  200  400  600  800 1000
## 
## $counts
## [1] 50  2  0  0  1
## 
## $density
## [1] 4.716981e-03 1.886792e-04 0.000000e+00 0.000000e+00 9.433962e-05
## 
## $mids
## [1] 100 300 500 700 900
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Alk
## $breaks
## [1]    0  500 1000 1500 2000 2500 3000 3500 4000
## 
## $counts
## [1] 13 25  7  1  6  0  0  1
## 
## $density
## [1] 4.905660e-04 9.433962e-04 2.641509e-04 3.773585e-05 2.264151e-04
## [6] 0.000000e+00 0.000000e+00 3.773585e-05
## 
## $mids
## [1]  250  750 1250 1750 2250 2750 3250 3750
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $TSS
## $breaks
## [1]   0  50 100 150 200 250
## 
## $counts
## [1] 31 18  3  0  1
## 
## $density
## [1] 0.0116981132 0.0067924528 0.0011320755 0.0000000000 0.0003773585
## 
## $mids
## [1]  25  75 125 175 225
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $PO4
## $breaks
## [1] 0.0e+00 5.0e-07 1.0e-06 1.5e-06 2.0e-06 2.5e-06 3.0e-06
## 
## $counts
## [1]  8 15  4  6 18  2
## 
## $density
## [1] 301886.8 566037.7 150943.4 226415.1 679245.3  75471.7
## 
## $mids
## [1] 2.50e-07 7.50e-07 1.25e-06 1.75e-06 2.25e-06 2.75e-06
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $NH4
## $breaks
##  [1] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040
## [10] 0.00045
## 
## $counts
## [1] 51  1  0  0  0  0  0  0  1
## 
## $density
## [1] 19245.2830   377.3585     0.0000     0.0000     0.0000     0.0000     0.0000
## [8]     0.0000   377.3585
## 
## $mids
## [1] 0.000025 0.000075 0.000125 0.000175 0.000225 0.000275 0.000325 0.000375
## [9] 0.000425
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $NO2
## $breaks
## [1] 0.0e+00 5.0e-06 1.0e-05 1.5e-05 2.0e-05
## 
## $counts
## [1] 49  2  1  1
## 
## $density
## [1] 184905.660   7547.170   3773.585   3773.585
## 
## $mids
## [1] 2.50e-06 7.50e-06 1.25e-05 1.75e-05
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $NO3
## $breaks
##  [1] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040
## [10] 0.00045
## 
## $counts
## [1] 21  6  9  6  4  2  2  1  2
## 
## $density
## [1] 7924.5283 2264.1509 3396.2264 2264.1509 1509.4340  754.7170  754.7170
## [8]  377.3585  754.7170
## 
## $mids
## [1] 0.000025 0.000075 0.000125 0.000175 0.000225 0.000275 0.000325 0.000375
## [9] 0.000425
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Br
## $breaks
## [1] 0.0e+00 5.0e-07 1.0e-06 1.5e-06 2.0e-06 2.5e-06
## 
## $counts
## [1] 18 22 10  0  3
## 
## $density
## [1] 679245.3 830188.7 377358.5      0.0 113207.5
## 
## $mids
## [1] 2.50e-07 7.50e-07 1.25e-06 1.75e-06 2.25e-06
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Cl
## $breaks
## [1] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040
## 
## $counts
## [1]  5 19 17  4  6  1  0  1
## 
## $density
## [1] 1886.7925 7169.8113 6415.0943 1509.4340 2264.1509  377.3585    0.0000
## [8]  377.3585
## 
## $mids
## [1] 0.000025 0.000075 0.000125 0.000175 0.000225 0.000275 0.000325 0.000375
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Fl
## $breaks
## [1] 1e-05 2e-05 3e-05 4e-05 5e-05 6e-05 7e-05
## 
## $counts
## [1] 10 28  8  4  1  2
## 
## $density
## [1] 18867.925 52830.189 15094.340  7547.170  1886.792  3773.585
## 
## $mids
## [1] 1.5e-05 2.5e-05 3.5e-05 4.5e-05 5.5e-05 6.5e-05
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $SO4
## $breaks
##  [1] 0e+00 1e-05 2e-05 3e-05 4e-05 5e-05 6e-05 7e-05 8e-05 9e-05
## 
## $counts
## [1]  3 34  9  3  1  2  0  0  1
## 
## $density
## [1]  5660.377 64150.943 16981.132  5660.377  1886.792  3773.585     0.000
## [8]     0.000  1886.792
## 
## $mids
## [1] 5.0e-06 1.5e-05 2.5e-05 3.5e-05 4.5e-05 5.5e-05 6.5e-05 7.5e-05 8.5e-05
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Na
## $breaks
## [1] 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012
## 
## $counts
## [1]  2 23 17  5  4  2
## 
## $density
## [1]  188.6792 2169.8113 1603.7736  471.6981  377.3585  188.6792
## 
## $mids
## [1] 0.0001 0.0003 0.0005 0.0007 0.0009 0.0011
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $K
## $breaks
## [1] 0e+00 2e-04 4e-04 6e-04 8e-04 1e-03
## 
## $counts
## [1] 39 12  1  0  1
## 
## $density
## [1] 3679.24528 1132.07547   94.33962    0.00000   94.33962
## 
## $mids
## [1] 1e-04 3e-04 5e-04 7e-04 9e-04
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Ca
## $breaks
##  [1] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040
## [10] 0.00045 0.00050 0.00055
## 
## $counts
##  [1]  2 23 17  6  0  0  3  1  0  0  1
## 
## $density
##  [1]  754.7170 8679.2453 6415.0943 2264.1509    0.0000    0.0000 1132.0755
##  [8]  377.3585    0.0000    0.0000  377.3585
## 
## $mids
##  [1] 0.000025 0.000075 0.000125 0.000175 0.000225 0.000275 0.000325 0.000375
##  [9] 0.000425 0.000475 0.000525
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Mg
## $breaks
## [1] 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030
## 
## $counts
## [1] 28 20  2  2  0  1
## 
## $density
## [1] 10566.0377  7547.1698   754.7170   754.7170     0.0000   377.3585
## 
## $mids
## [1] 0.000025 0.000075 0.000125 0.000175 0.000225 0.000275
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $TDN
## $breaks
## [1] 0e+00 1e-04 2e-04 3e-04 4e-04 5e-04 6e-04
## 
## $counts
## [1] 26 11  9  4  2  1
## 
## $density
## [1] 4905.6604 2075.4717 1698.1132  754.7170  377.3585  188.6792
## 
## $mids
## [1] 0.00005 0.00015 0.00025 0.00035 0.00045 0.00055
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
wc[,-2]<-log(wc[,-2]) #log-transform concentration data (not pH)
apply(wc, 2, hist) #check normal distribution

## $cond
## $breaks
## [1] 3.5 4.0 4.5 5.0 5.5 6.0 6.5
## 
## $counts
## [1]  4 20 17 11  0  1
## 
## $density
## [1] 0.15094340 0.75471698 0.64150943 0.41509434 0.00000000 0.03773585
## 
## $mids
## [1] 3.75 4.25 4.75 5.25 5.75 6.25
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $pH
## $breaks
## [1] 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8
## 
## $counts
## [1]  3 11 18 13  6  0  1  1
## 
## $density
## [1] 0.28301887 1.03773585 1.69811321 1.22641509 0.56603774 0.00000000 0.09433962
## [8] 0.09433962
## 
## $mids
## [1] 7.3 7.5 7.7 7.9 8.1 8.3 8.5 8.7
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $turb
## $breaks
## [1] -4 -2  0  2  4  6  8
## 
## $counts
## [1] 21  1  3 14 13  1
## 
## $density
## [1] 0.198113208 0.009433962 0.028301887 0.132075472 0.122641509 0.009433962
## 
## $mids
## [1] -3 -1  1  3  5  7
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Alk
## $breaks
## [1] 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5
## 
## $counts
## [1]  1  1 24 14  6  6  1
## 
## $density
## [1] 0.03773585 0.03773585 0.90566038 0.52830189 0.22641509 0.22641509 0.03773585
## 
## $mids
## [1] 5.25 5.75 6.25 6.75 7.25 7.75 8.25
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $TSS
## $breaks
## [1] 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5
## 
## $counts
## [1]  1  6 11  7 10 14  3  1
## 
## $density
## [1] 0.03773585 0.22641509 0.41509434 0.26415094 0.37735849 0.52830189 0.11320755
## [8] 0.03773585
## 
## $mids
## [1] 1.75 2.25 2.75 3.25 3.75 4.25 4.75 5.25
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $PO4
## $breaks
## [1] -16.0 -15.5 -15.0 -14.5 -14.0 -13.5 -13.0 -12.5
## 
## $counts
## [1]  5  1  2 10  8 18  9
## 
## $density
## [1] 0.18867925 0.03773585 0.07547170 0.37735849 0.30188679 0.67924528 0.33962264
## 
## $mids
## [1] -15.75 -15.25 -14.75 -14.25 -13.75 -13.25 -12.75
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $NH4
## $breaks
## [1] -15 -14 -13 -12 -11 -10  -9  -8  -7
## 
## $counts
## [1]  2  0  4 19 25  2  0  1
## 
## $density
## [1] 0.03773585 0.00000000 0.07547170 0.35849057 0.47169811 0.03773585 0.00000000
## [8] 0.01886792
## 
## $mids
## [1] -14.5 -13.5 -12.5 -11.5 -10.5  -9.5  -8.5  -7.5
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $NO2
## $breaks
## [1] -16 -15 -14 -13 -12 -11 -10
## 
## $counts
## [1]  4 17 24  4  3  1
## 
## $density
## [1] 0.07547170 0.32075472 0.45283019 0.07547170 0.05660377 0.01886792
## 
## $mids
## [1] -15.5 -14.5 -13.5 -12.5 -11.5 -10.5
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $NO3
## $breaks
##  [1] -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0  -9.5  -9.0  -8.5  -8.0  -7.5
## 
## $counts
##  [1]  1  0  2  0 11  6  6  2 15  7  3
## 
## $density
##  [1] 0.03773585 0.00000000 0.07547170 0.00000000 0.41509434 0.22641509
##  [7] 0.22641509 0.07547170 0.56603774 0.26415094 0.11320755
## 
## $mids
##  [1] -12.75 -12.25 -11.75 -11.25 -10.75 -10.25  -9.75  -9.25  -8.75  -8.25
## [11]  -7.75
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Br
## $breaks
## [1] -16.0 -15.5 -15.0 -14.5 -14.0 -13.5 -13.0 -12.5
## 
## $counts
## [1]  2  6 10 16 14  3  2
## 
## $density
## [1] 0.0754717 0.2264151 0.3773585 0.6037736 0.5283019 0.1132075 0.0754717
## 
## $mids
## [1] -15.75 -15.25 -14.75 -14.25 -13.75 -13.25 -12.75
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Cl
## $breaks
##  [1] -12.5 -12.0 -11.5 -11.0 -10.5 -10.0  -9.5  -9.0  -8.5  -8.0  -7.5
## 
## $counts
##  [1]  1  0  0  0  4  8 23 10  6  1
## 
## $density
##  [1] 0.03773585 0.00000000 0.00000000 0.00000000 0.15094340 0.30188679
##  [7] 0.86792453 0.37735849 0.22641509 0.03773585
## 
## $mids
##  [1] -12.25 -11.75 -11.25 -10.75 -10.25  -9.75  -9.25  -8.75  -8.25  -7.75
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Fl
## $breaks
##  [1] -11.6 -11.4 -11.2 -11.0 -10.8 -10.6 -10.4 -10.2 -10.0  -9.8  -9.6
## 
## $counts
##  [1]  1  2  2  8 12 14  3  7  2  2
## 
## $density
##  [1] 0.09433962 0.18867925 0.18867925 0.75471698 1.13207547 1.32075472
##  [7] 0.28301887 0.66037736 0.18867925 0.18867925
## 
## $mids
##  [1] -11.5 -11.3 -11.1 -10.9 -10.7 -10.5 -10.3 -10.1  -9.9  -9.7
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $SO4
## $breaks
## [1] -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0  -9.5  -9.0
## 
## $counts
## [1]  3  0  0 21 22  3  3  1
## 
## $density
## [1] 0.11320755 0.00000000 0.00000000 0.79245283 0.83018868 0.11320755 0.11320755
## [8] 0.03773585
## 
## $mids
## [1] -12.75 -12.25 -11.75 -11.25 -10.75 -10.25  -9.75  -9.25
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Na
## $breaks
## [1] -10.0  -9.5  -9.0  -8.5  -8.0  -7.5  -7.0  -6.5
## 
## $counts
## [1]  1  0  1 17 21  9  4
## 
## $density
## [1] 0.03773585 0.00000000 0.03773585 0.64150943 0.79245283 0.33962264 0.15094340
## 
## $mids
## [1] -9.75 -9.25 -8.75 -8.25 -7.75 -7.25 -6.75
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $K
## $breaks
## [1] -12 -11 -10  -9  -8  -7  -6
## 
## $counts
## [1]  1  0 12 36  3  1
## 
## $density
## [1] 0.01886792 0.00000000 0.22641509 0.67924528 0.05660377 0.01886792
## 
## $mids
## [1] -11.5 -10.5  -9.5  -8.5  -7.5  -6.5
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Ca
## $breaks
## [1] -11.0 -10.5 -10.0  -9.5  -9.0  -8.5  -8.0  -7.5
## 
## $counts
## [1]  1  0 12 23 12  3  2
## 
## $density
## [1] 0.03773585 0.00000000 0.45283019 0.86792453 0.45283019 0.11320755 0.07547170
## 
## $mids
## [1] -10.75 -10.25  -9.75  -9.25  -8.75  -8.25  -7.75
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $Mg
## $breaks
## [1] -12.0 -11.5 -11.0 -10.5 -10.0  -9.5  -9.0  -8.5  -8.0
## 
## $counts
## [1]  1  0  1 21 22  3  4  1
## 
## $density
## [1] 0.03773585 0.00000000 0.03773585 0.79245283 0.83018868 0.11320755 0.15094340
## [8] 0.03773585
## 
## $mids
## [1] -11.75 -11.25 -10.75 -10.25  -9.75  -9.25  -8.75  -8.25
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"
## 
## $TDN
## $breaks
## [1] -11.0 -10.5 -10.0  -9.5  -9.0  -8.5  -8.0  -7.5  -7.0
## 
## $counts
## [1]  3 13  6  5 10 11  4  1
## 
## $density
## [1] 0.11320755 0.49056604 0.22641509 0.18867925 0.37735849 0.41509434 0.15094340
## [8] 0.03773585
## 
## $mids
## [1] -10.75 -10.25  -9.75  -9.25  -8.75  -8.25  -7.75  -7.25
## 
## $xname
## [1] "newX[, i]"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

plot(wc) # check potential for PCA by correlation plot


pca<-prcomp(scale(wc),center=F,scale.=F)
pca<-prcomp(wc,center=T,scale.=T)   # equivalent to line above

summary(pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.8392 1.8644 1.22808 1.10689 0.84577 0.80664 0.78730
## Proportion of Variance 0.4478 0.1931 0.08379 0.06807 0.03974 0.03615 0.03444
## Cumulative Proportion  0.4478 0.6410 0.72474 0.79281 0.83255 0.86870 0.90314
##                            PC8    PC9    PC10    PC11   PC12    PC13    PC14
## Standard deviation     0.66233 0.6392 0.53771 0.40863 0.3445 0.32772 0.29622
## Proportion of Variance 0.02437 0.0227 0.01606 0.00928 0.0066 0.00597 0.00487
## Cumulative Proportion  0.92751 0.9502 0.96627 0.97554 0.9821 0.98810 0.99298
##                           PC15    PC16   PC17    PC18
## Standard deviation     0.25252 0.21425 0.1038 0.07711
## Proportion of Variance 0.00354 0.00255 0.0006 0.00033
## Cumulative Proportion  0.99652 0.99907 0.9997 1.00000

pca$sdev # stdevs of PCs (squares are eigenvalues); eigenvalues correspond to the variance of their respective principal component (PC)
##  [1] 2.83921334 1.86442507 1.22808376 1.10689311 0.84576999 0.80663510
##  [7] 0.78729929 0.66233089 0.63915683 0.53770934 0.40862593 0.34454651
## [13] 0.32772103 0.29622209 0.25251770 0.21424557 0.10379779 0.07711088
eigenvals(pca)==pca$sdev^2
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 
## PC17 PC18 
## TRUE TRUE

head(scores<-pca$x) # site scores on all PCs
##         PC1      PC2         PC3        PC4        PC5        PC6         PC7
## 1 -4.749326 3.557587  0.52830345  1.6885404  0.4637009  0.6815660  0.06534355
## 2 -3.830547 3.155164  0.07572526 -2.0332541  0.2644377  1.1129832 -0.55448962
## 3 -3.267373 4.108748 -0.32897368  0.7711051 -1.0759520 -0.2527134  0.44879703
## 4 -3.680797 2.292490  0.80951314  0.2051344 -0.2780356 -0.6345223  0.05431164
## 5 -8.999719 1.602003 -3.03923365 -1.2132318 -1.0636252 -1.2871947  0.42462524
## 7  2.315295 1.447056  0.53769126  0.3896117  0.5148749 -1.6163461 -0.06189673
##          PC8        PC9        PC10       PC11         PC12        PC13
## 1  0.9056796 -0.1326339 -0.36156648  0.3464919 -0.343851835 -0.09665447
## 2  1.2660092 -0.0486270 -1.04265251 -0.2685484 -0.237503649 -0.68256243
## 3 -0.1695538 -0.4501332 -0.45996058  0.5662810  0.124965534 -0.53160488
## 4 -0.5018929 -0.9758938 -0.49191459  0.1223105 -0.230524069  0.50704177
## 5  0.3163298 -0.1965387 -0.05913828 -0.1130011  1.000404068  0.11682215
## 7 -1.1804198 -1.0651942 -0.46082573 -1.2090894 -0.008212293 -0.32546687
##          PC14       PC15         PC16        PC17        PC18
## 1  0.36104029  0.1178840  0.004293307  0.04989172  0.05965282
## 2 -0.14602705 -0.1922231  0.155271851  0.28350624 -0.09569322
## 3 -0.40743689  0.4087700  0.288081091 -0.24071374  0.04519902
## 4 -0.01593147 -0.3959628 -0.098755194 -0.01316164  0.08656537
## 5 -0.17548689  0.1746190 -0.311008889  0.19590434 -0.02083158
## 7  0.58480671 -0.4687339 -0.016056662 -0.03962186  0.01500699

head(scale(wc) %*% pca$rotation) # to manually compute scores from variables and loadings
##         PC1      PC2         PC3        PC4        PC5        PC6         PC7
## 1 -4.749326 3.557587  0.52830345  1.6885404  0.4637009  0.6815660  0.06534355
## 2 -3.830547 3.155164  0.07572526 -2.0332541  0.2644377  1.1129832 -0.55448962
## 3 -3.267373 4.108748 -0.32897368  0.7711051 -1.0759520 -0.2527134  0.44879703
## 4 -3.680797 2.292490  0.80951314  0.2051344 -0.2780356 -0.6345223  0.05431164
## 5 -8.999719 1.602003 -3.03923365 -1.2132318 -1.0636252 -1.2871947  0.42462524
## 7  2.315295 1.447056  0.53769126  0.3896117  0.5148749 -1.6163461 -0.06189673
##          PC8        PC9        PC10       PC11         PC12        PC13
## 1  0.9056796 -0.1326339 -0.36156648  0.3464919 -0.343851835 -0.09665447
## 2  1.2660092 -0.0486270 -1.04265251 -0.2685484 -0.237503649 -0.68256243
## 3 -0.1695538 -0.4501332 -0.45996058  0.5662810  0.124965534 -0.53160488
## 4 -0.5018929 -0.9758938 -0.49191459  0.1223105 -0.230524069  0.50704177
## 5  0.3163298 -0.1965387 -0.05913828 -0.1130011  1.000404068  0.11682215
## 7 -1.1804198 -1.0651942 -0.46082573 -1.2090894 -0.008212293 -0.32546687
##          PC14       PC15         PC16        PC17        PC18
## 1  0.36104029  0.1178840  0.004293307  0.04989172  0.05965282
## 2 -0.14602705 -0.1922231  0.155271851  0.28350624 -0.09569322
## 3 -0.40743689  0.4087700  0.288081091 -0.24071374  0.04519902
## 4 -0.01593147 -0.3959628 -0.098755194 -0.01316164  0.08656537
## 5 -0.17548689  0.1746190 -0.311008889  0.19590434 -0.02083158
## 7  0.58480671 -0.4687339 -0.016056662 -0.03962186  0.01500699
head(loadings<-pca$rotation) # the variable loadings
##              PC1        PC2        PC3         PC4         PC5         PC6
## cond -0.26935347  0.1438792 -0.2841620 -0.03056776  0.12191093  0.08071683
## pH    0.04777530  0.1215468  0.3562802 -0.59539396 -0.49268822 -0.29216145
## turb -0.14155181 -0.3894707 -0.2177088 -0.32746528 -0.03611414  0.24629283
## Alk  -0.17098754  0.3503345 -0.1962747  0.09297582  0.10935101 -0.12209184
## TSS  -0.14508786 -0.3207217 -0.2732119 -0.44573924 -0.03791156  0.26821922
## PO4  -0.01446834 -0.2443586  0.4595895 -0.23241378  0.72222431 -0.10638692
##              PC7         PC8         PC9         PC10        PC11         PC12
## cond -0.41973560 -0.17562478 -0.09881877 -0.006566902  0.64824153 -0.064479508
## pH   -0.22576487 -0.26192817  0.06281913 -0.192775926  0.05253252 -0.003739136
## turb -0.04788218  0.07002402 -0.01614397  0.011308932 -0.32201794 -0.137669971
## Alk  -0.42805174 -0.15744036 -0.37748332 -0.112749767 -0.59639437  0.187029632
## TSS   0.02529992  0.16003510 -0.08747407 -0.150755304 -0.03128159  0.262605089
## PO4  -0.05401838 -0.17326809 -0.05875984 -0.015299767 -0.07294578 -0.177318022
##             PC13        PC14        PC15         PC16          PC17
## cond -0.12739932  0.02559159 -0.11856620 -0.336469042 -0.0003470607
## pH   -0.04129351  0.07363790  0.02360840  0.043866907 -0.0064374198
## turb -0.59531621 -0.04723956  0.24918988 -0.223721709 -0.1124305273
## Alk   0.01598820  0.05220398 -0.07854688  0.058152265  0.0221005482
## TSS   0.55485685 -0.08756056 -0.29311696  0.002605816  0.0204648859
## PO4   0.03082657  0.13171082 -0.19384168 -0.090576799  0.0664719576
##              PC18
## cond  0.135379479
## pH    0.017462207
## turb  0.005186753
## Alk   0.084188983
## TSS   0.030558062
## PO4  -0.005485131


# how many axes should be kept?
#there is not the one well-accepted method
#1. visual examination
#looking for a point at which the proportion of variance explained by each subsequent principal component drops off (elbow of scree plot)
screeplot(pca, npcs=length(pca$sdev), type = "lines") # plots eigenvalues vs. principal components (PCs) #
#2 Kaiser-Guttman criterion
# retain PCs with eigenvalues>1, PCs with eigenvalues<1 contain less informatin than one original variable
abline(h=1,col="red")

# eigenvalue of PC5<1
#3. variance explained method: set an "arbitrary" threshold of explained variance (e.g. 70%, 80%, 90%) and choose as many PCs till threshold is reached.
summary(pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.8392 1.8644 1.22808 1.10689 0.84577 0.80664 0.78730
## Proportion of Variance 0.4478 0.1931 0.08379 0.06807 0.03974 0.03615 0.03444
## Cumulative Proportion  0.4478 0.6410 0.72474 0.79281 0.83255 0.86870 0.90314
##                            PC8    PC9    PC10    PC11   PC12    PC13    PC14
## Standard deviation     0.66233 0.6392 0.53771 0.40863 0.3445 0.32772 0.29622
## Proportion of Variance 0.02437 0.0227 0.01606 0.00928 0.0066 0.00597 0.00487
## Cumulative Proportion  0.92751 0.9502 0.96627 0.97554 0.9821 0.98810 0.99298
##                           PC15    PC16   PC17    PC18
## Standard deviation     0.25252 0.21425 0.1038 0.07711
## Proportion of Variance 0.00354 0.00255 0.0006 0.00033
## Cumulative Proportion  0.99652 0.99907 0.9997 1.00000

# -> the first 3-4 PCs seem useful, and just PC1 and PC2 alone are already explaining a lot of overall variance


# for a DISTANCE BIPLOT (focus is on sites, "scaling 1")
# each principal component has variance given by eigenvalue, loadings remain unscaled
str(landuse)
##  chr [1:53] "A" "A" "A" "A" "A" "M" "F" "M" "M" "M" "F" "A" "A" "A" "A" "A" ...
landuse<-as.factor(landuse)
plot(scores[,1:2],asp=1,pch=21, cex=2, bg=landuse) # note asp=1
arrows<-loadings*15 # with extension factor
Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="darkgreen")
text(x=arrows[,1]*1.2,y=arrows[,2]*1.2,labels=names(wc),cex=0.8)
#text(scores[,1:2], rownames(wc), pos=2)
legend("bottomright", legend=c("Agriculture", "Mixed", "Forest"), pch=21, pt.cex=2, pt.bg=c("black","green","darkred"), cex=1.5)


biplot(pca,scale=0)

# in this plot:
# 1) distances among sites are approximating true Euclidean distances in multivariate space
# 2) angles between arrows do not reflect correlations among variables
# 3) projecting site on descriptor at right angle gives its appr. descriptor value

# for a CORRELATION BIPLOT (focus is on variables, "scaling 2")
# each principal component is weighted by 1/sqrt(eigenvalue), so it has variance 1
var(scores[,1]/pca$sdev[1]) # just demo
## [1] 1
plot(scores[,1]/pca$sdev[1],scores[,2]/pca$sdev[2],pch=21,bg=landuse,asp=1)
# loadings are weighted by sqrt(eigenvalues)
arrows<-loadings*matrix(pca$sdev,nrow=nrow(loadings),ncol=ncol(loadings),byrow=TRUE)
arrows<-arrows*2 # choose extension factor
Arrows(x0=0,y0=0,x1=arrows[,1],y1=arrows[,2],col="purple")
# as alternative just compute correlation of scores with original data ("structure coefficients")
(structure<-cor(wc,scores))
##              PC1         PC2         PC3         PC4          PC5          PC6
## cond -0.76475195  0.26825195 -0.34897469 -0.03383524  0.103108608  0.065109031
## pH    0.13564426  0.22661496  0.43754197 -0.65903747 -0.416700907 -0.235667684
## turb -0.40189578 -0.72613887 -0.26736468 -0.36246906 -0.030544259  0.198668445
## Alk  -0.48547011  0.65317251 -0.24104174  0.10291430  0.092485802 -0.098483565
## TSS  -0.41193538 -0.59796155 -0.33552705 -0.49338569 -0.032064459  0.216355039
## PO4  -0.04107869 -0.45558830  0.56441437 -0.25725721  0.610835646 -0.085815422
## NH4  -0.59306916  0.02901967 -0.52211239 -0.05077205  0.146005150 -0.407246534
## NO2  -0.60159868 -0.44979709 -0.11708309 -0.06542819 -0.083433731 -0.478631077
## NO3  -0.51894312 -0.68971906  0.17118492  0.35407824 -0.145546662 -0.008334060
## Br   -0.92514493 -0.10433460  0.06301001  0.05315603  0.008796639  0.221187748
## Cl   -0.82983473 -0.12663016  0.38234130  0.16849740  0.028517422 -0.093976441
## Fl   -0.85617616 -0.08722719  0.14001888 -0.03859390  0.108893567 -0.027198849
## SO4  -0.81130553  0.18629947  0.30840708  0.17786031 -0.164712321  0.009860507
## Na   -0.88109521  0.13206079  0.10468309  0.04519175 -0.161228683  0.142295300
## K    -0.79637284  0.37181544  0.19278724 -0.11851266  0.026616893  0.075563024
## Ca   -0.73330911  0.56721172 -0.04684159 -0.17535185  0.084387922  0.046307940
## Mg   -0.76241145  0.49212149  0.03837621 -0.19104400  0.036005859  0.052768803
## TDN  -0.66409999 -0.60500271 -0.04421175  0.24134387 -0.151522130 -0.055259365
##              PC7         PC8          PC9         PC10         PC11
## cond -0.33045754 -0.11632172 -0.063160692 -0.003531084  0.264888300
## pH   -0.17774452 -0.17348312  0.040151278 -0.103657416  0.021466149
## turb -0.03769761  0.04637907 -0.010318528  0.006080919 -0.131584882
## Alk  -0.33700483 -0.10427761 -0.241271045 -0.060626603 -0.243702206
## TSS   0.01991861  0.10599619 -0.055909646 -0.081062536 -0.012782470
## PO4  -0.04252863 -0.11476081 -0.037556751 -0.008226828 -0.029807536
## NH4   0.14949569 -0.04291440  0.365863964 -0.155284000 -0.033708442
## NO2   0.07756756  0.19275231 -0.298888832  0.209087387  0.052794089
## NO3   0.01371875 -0.25604120 -0.048509024 -0.033794332 -0.013701263
## Br    0.04173239  0.00127664 -0.052263599 -0.162136464  0.069169573
## Cl   -0.03699389  0.13621525  0.008198462 -0.180611063  0.032692795
## Fl   -0.33251999  0.16703916  0.191260615  0.124728010  0.003445899
## SO4   0.05600956  0.30157467 -0.015244843 -0.159709162 -0.046927081
## Na   -0.13133773 -0.01459469  0.187683610  0.211586965 -0.073979848
## K     0.21133725 -0.08141740  0.106111667  0.223653004 -0.011006163
## Ca    0.21377568 -0.09313882 -0.149458288 -0.054964294  0.017660115
## Mg    0.32651015 -0.06474423 -0.077876661  0.034731711 -0.011551626
## TDN   0.04747067 -0.30977212 -0.004718641  0.007756368 -0.019322840
##              PC12         PC13          PC14         PC15          PC16
## cond -0.022216190 -0.041751438  0.0075807958 -0.029940066 -0.0720870031
## pH   -0.001288306 -0.013532752  0.0218131727  0.005961538  0.0093982907
## turb -0.047433709 -0.195097641 -0.0139934007  0.062924857 -0.0479313860
## Alk   0.064440408  0.005239670  0.0154639711 -0.019834478  0.0124588654
## TSS   0.090479668  0.181838259 -0.0259373713 -0.074017223  0.0005582845
## PO4  -0.061094306  0.010102514  0.0390156553 -0.048948456 -0.0194056782
## NH4  -0.019558805 -0.008153826  0.0261212836 -0.022615617  0.0068037616
## NO2  -0.018137139 -0.022027853  0.0217333801 -0.022298082  0.0303642603
## NO3  -0.015197048  0.066748065  0.0090589510  0.018038656 -0.0081950848
## Br    0.003987132 -0.074486267  0.1268513291 -0.005079333  0.1418328916
## Cl    0.108916049 -0.077892719 -0.1965043400 -0.029094888  0.0066897381
## Fl    0.021179746  0.098884906  0.0044583717  0.154701308  0.0292928111
## SO4  -0.064893733  0.034208880  0.1163038494 -0.018049823 -0.1085844000
## Na   -0.153097268  0.003712544 -0.0563765190 -0.126996360  0.0355530451
## K     0.216803483 -0.052579706  0.0705418244 -0.028508949 -0.0372784074
## Ca   -0.081684815  0.049612428 -0.0519750553  0.054216175 -0.0041896038
## Mg   -0.068998656  0.024034873 -0.0651154780  0.058400432  0.0043412431
## TDN   0.023214045  0.028784640 -0.0001785383  0.029159768 -0.0353732943
##               PC17          PC18
## cond -3.602414e-05  0.0104392312
## pH   -6.681899e-04  0.0013465262
## turb -1.167004e-02  0.0003999551
## Alk   2.293988e-03  0.0064918868
## TSS   2.124210e-03  0.0023563591
## PO4   6.899642e-03 -0.0004229633
## NH4  -9.032968e-03  0.0036488485
## NO2  -2.342515e-03 -0.0016575665
## NO3  -6.145284e-02  0.0262899055
## Br    8.706134e-03 -0.0017381812
## Cl    6.151965e-04 -0.0017187344
## Fl    6.209194e-04 -0.0016118828
## SO4   6.993261e-03 -0.0001298039
## Na   -2.090597e-04 -0.0059058664
## K    -1.630078e-02 -0.0026053297
## Ca   -3.672943e-02 -0.0461913573
## Mg    3.246239e-02  0.0471881766
## TDN   6.262864e-02 -0.0258731171
structure<-2*structure
Arrows(x0=0,y0=0,x1=structure[,1],y1=structure[,2],col="red")
text(x=arrows[,1]*1.3,y=arrows[,2]*1.2,labels=names(wc),cex=0.7)


biplot(pca,scale=1)

# in this plot
# 1) distances among sites are not approximating true Euclidean distances in multivariate space
# 2) angles between arrows reflect correlations among variables (NOT proximity of arrow heads)
# 3) projecting site on descriptor at right angle gives its appr. descriptor value

# PCA using the rda function from the vegan package
pca2<-rda(X=wc,scale=TRUE)
summary(pca2,scaling=1)
## 
## Call:
## rda(X = wc, scale = TRUE) 
## 
## Partitioning of correlations:
##               Inertia Proportion
## Total              18          1
## Unconstrained      18          1
## 
## Eigenvalues, and their contribution to the correlations 
## 
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Eigenvalue            8.0611 3.4761 1.50819 1.22521 0.71533 0.65066 0.61984
## Proportion Explained  0.4478 0.1931 0.08379 0.06807 0.03974 0.03615 0.03444
## Cumulative Proportion 0.4478 0.6410 0.72474 0.79281 0.83255 0.86870 0.90314
##                           PC8    PC9    PC10     PC11     PC12     PC13
## Eigenvalue            0.43868 0.4085 0.28913 0.166975 0.118712 0.107401
## Proportion Explained  0.02437 0.0227 0.01606 0.009276 0.006595 0.005967
## Cumulative Proportion 0.92751 0.9502 0.96627 0.975542 0.982137 0.988104
##                           PC14     PC15    PC16      PC17      PC18
## Eigenvalue            0.087748 0.063765 0.04590 0.0107740 0.0059461
## Proportion Explained  0.004875 0.003543 0.00255 0.0005986 0.0003303
## Cumulative Proportion 0.992979 0.996521 0.99907 0.9996697 1.0000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  5.531195 
## 
## 
## Species scores
## 
##           PC1      PC2     PC3     PC4      PC5      PC6
## cond -1.48985  0.79582 -1.5718 -0.1691  0.67431  0.44646
## pH    0.26425  0.67230  1.9707 -3.2932 -2.72515 -1.61600
## turb -0.78295 -2.15424 -1.2042 -1.8113 -0.19975  1.36229
## Alk  -0.94577  1.93777 -1.0856  0.5143  0.60484 -0.67531
## TSS  -0.80251 -1.77397 -1.5112 -2.4655 -0.20970  1.48357
## PO4  -0.08003 -1.35160  2.5421 -1.2855  3.99476 -0.58845
## NH4  -1.15538  0.08609 -2.3516 -0.2537  0.95485 -2.79254
## NO2  -1.17200 -1.33441 -0.5273 -0.3269 -0.54564 -3.28203
## NO3  -1.01098 -2.04619  0.7710  1.7693 -0.95185 -0.05715
## Br   -1.80232 -0.30953  0.2838  0.2656  0.05753  1.51671
## Cl   -1.61664 -0.37567  1.7220  0.8420  0.18650 -0.64441
## Fl   -1.66795 -0.25878  0.6306 -0.1929  0.71215 -0.18651
## SO4  -1.58054  0.55270  1.3890  0.8888 -1.07719  0.06761
## Na   -1.71650  0.39179  0.4715  0.2258 -1.05441  0.97574
## K    -1.55145  1.10307  0.8683 -0.5922  0.17407  0.51814
## Ca   -1.42859  1.68275 -0.2110 -0.8762  0.55188  0.31754
## Mg   -1.48529  1.45998  0.1728 -0.9547  0.23547  0.36184
## TDN  -1.29376 -1.79486 -0.1991  1.2060 -0.99093 -0.37892
## 
## 
## Site scores (weighted sums of species scores)
## 
##          PC1      PC2       PC3      PC4       PC5       PC6
## 1  -0.858644  0.64319  0.095513  0.30528  0.083834  0.123222
## 2  -0.692535  0.57043  0.013691 -0.36760  0.047808  0.201219
## 3  -0.590717  0.74283 -0.059476  0.13941 -0.194524 -0.045689
## 4  -0.665461  0.41447  0.146354  0.03709 -0.050267 -0.114717
## 5  -1.627084  0.28963 -0.549471 -0.21934 -0.192296 -0.232715
## 7   0.418588  0.26162  0.097211  0.07044  0.093086 -0.292224
## 8   0.990156  0.31409 -0.773129 -0.34073 -0.067209  0.150470
## 9  -0.298376  0.20548 -0.191656  0.70420 -0.194790  0.113638
## 10 -0.509697  0.18895  0.022804 -0.32743  0.211013  0.044924
## 11  0.158309  0.31689 -0.173319  0.16996  0.237407  0.066935
## 12  0.764213  0.20918 -0.027062  0.06540 -0.147885 -0.096221
## 13 -0.178540 -0.35882 -0.150284 -0.02218  0.196198  0.030318
## 14  1.242212 -0.56161 -0.617922  0.08199  0.106896 -0.193683
## 15 -0.368970 -0.58101 -0.007750 -0.05095 -0.062378  0.057150
## 16 -0.669563 -0.04282 -0.047768  0.21397  0.152919 -0.277781
## 17 -0.761438 -0.24526  0.029662  0.03186  0.056651 -0.066478
## 18 -0.425436 -0.04832 -0.024151  0.02310  0.017175  0.226641
## 19 -0.234255 -0.27163  0.068958 -0.19446 -0.054914  0.018765
## 20  0.028982 -0.19608 -0.178394  0.05594 -0.199762  0.196180
## 21 -0.338971 -0.32675  0.150791 -0.08567 -0.105597 -0.022958
## 22 -0.128236 -0.32483  0.196796 -0.01866 -0.024355  0.385547
## 23 -0.363328 -0.18509  0.059418 -0.09107  0.017646 -0.048482
## 24 -0.596633 -0.45527  0.025255 -0.01668  0.055688 -0.127478
## 25 -0.030951 -0.36593  0.256660 -0.43481 -0.257818 -0.221803
## 26  0.192175 -0.44357  0.021972 -0.10034  0.052713  0.006413
## 27  0.041707 -0.25523 -0.121742  0.14506 -0.148207  0.038812
## 28  0.076520 -0.45121 -0.037076 -0.03006  0.096141 -0.035485
## 29  0.088275 -0.39706  0.026841  0.03825  0.115360  0.051904
## 30 -0.162040 -0.39082 -0.260650  0.12857 -0.041694  0.066871
## 31  0.400213  0.29276  0.183117  0.16617 -0.022077  0.210915
## 32  0.194940  0.27892  0.011061  0.10721  0.008114 -0.041581
## 33  0.096061 -0.04411  0.261267  0.44553 -0.062379 -0.335050
## 34 -0.082955  0.21771  0.144274 -0.37113  0.087125 -0.031603
## 35 -0.004696 -0.19546  0.173180 -0.07966  0.013836 -0.003374
## 36  0.217794  0.35991 -0.002031 -0.01349 -0.044245 -0.008067
## 37  0.530344  0.16473 -0.072752  0.08482 -0.179419  0.092896
## 38  0.386837  0.15010  0.163132 -0.04098  0.156511 -0.096270
## 39  0.454369  0.32809  0.158876  0.05545  0.229962  0.159861
## 40  0.370033  0.18849  0.150137 -0.01702  0.220776 -0.060391
## 41  0.733763  0.09891  0.543673 -0.25807 -0.215391 -0.102342
## 42  0.751397  0.28756 -0.018136  0.04335 -0.448897 -0.150765
## 43  0.260255  0.27344  0.283415 -0.15403  0.134286  0.053973
## 44  0.626689  0.03673  0.109380  0.08025  0.171751 -0.177150
## 45  0.422041  0.66969 -0.305212 -0.28421  0.194072 -0.040262
## 46 -0.308450 -0.38097 -0.156424 -0.03333  0.027181 -0.080304
## 47 -0.090387 -0.08771 -0.051645 -0.04502  0.049220 -0.105366
## 48  0.243109  0.17544  0.055587  0.16297  0.129623 -0.029377
## 49 -0.098095 -0.03125 -0.176481 -0.08163 -0.086694  0.113537
## 50 -0.123741 -0.28603 -0.028502  0.04567 -0.089356  0.096000
## 51 -0.220621 -0.11869  0.207751  0.06654 -0.025200  0.046210
## 52  0.447317 -0.04547  0.241645 -0.09827 -0.339366  0.198683
## 53  0.154271 -0.24656  0.021007  0.13312  0.174027  0.094545
## 54  0.139250 -0.34166  0.111607  0.17523  0.117702  0.191989
summary(pca2,scaling=2)
## 
## Call:
## rda(X = wc, scale = TRUE) 
## 
## Partitioning of correlations:
##               Inertia Proportion
## Total              18          1
## Unconstrained      18          1
## 
## Eigenvalues, and their contribution to the correlations 
## 
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Eigenvalue            8.0611 3.4761 1.50819 1.22521 0.71533 0.65066 0.61984
## Proportion Explained  0.4478 0.1931 0.08379 0.06807 0.03974 0.03615 0.03444
## Cumulative Proportion 0.4478 0.6410 0.72474 0.79281 0.83255 0.86870 0.90314
##                           PC8    PC9    PC10     PC11     PC12     PC13
## Eigenvalue            0.43868 0.4085 0.28913 0.166975 0.118712 0.107401
## Proportion Explained  0.02437 0.0227 0.01606 0.009276 0.006595 0.005967
## Cumulative Proportion 0.92751 0.9502 0.96627 0.975542 0.982137 0.988104
##                           PC14     PC15    PC16      PC17      PC18
## Eigenvalue            0.087748 0.063765 0.04590 0.0107740 0.0059461
## Proportion Explained  0.004875 0.003543 0.00255 0.0005986 0.0003303
## Cumulative Proportion 0.992979 0.996521 0.99907 0.9996697 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  5.531195 
## 
## 
## Species scores
## 
##           PC1      PC2      PC3      PC4      PC5      PC6
## cond -0.99702  0.34972 -0.45496 -0.04411  0.13442  0.08488
## pH    0.17684  0.29544  0.57043 -0.85920 -0.54326 -0.30724
## turb -0.52396 -0.94668 -0.34857 -0.47256 -0.03982  0.25901
## Alk  -0.63291  0.85155 -0.31425  0.13417  0.12058 -0.12839
## TSS  -0.53705 -0.77957 -0.43743 -0.64323 -0.04180  0.28207
## PO4  -0.05355 -0.59396  0.73584 -0.33539  0.79636 -0.11188
## NH4  -0.77319  0.03783 -0.68069 -0.06619  0.19035 -0.53093
## NO2  -0.78431 -0.58641 -0.15264 -0.08530 -0.10877 -0.62400
## NO3  -0.67655 -0.89920  0.22318  0.46162 -0.18975 -0.01087
## Br   -1.20613 -0.13602  0.08215  0.06930  0.01147  0.28837
## Cl   -1.08187 -0.16509  0.49846  0.21967  0.03718 -0.12252
## Fl   -1.11621 -0.11372  0.18254 -0.05032  0.14197 -0.03546
## SO4  -1.05771  0.24288  0.40207  0.23188 -0.21474  0.01286
## Na   -1.14870  0.17217  0.13648  0.05892 -0.21020  0.18551
## K    -1.03824  0.48474  0.25134 -0.15451  0.03470  0.09851
## Ca   -0.95603  0.73948 -0.06107 -0.22861  0.11002  0.06037
## Mg   -0.99397  0.64159  0.05003 -0.24907  0.04694  0.06880
## TDN  -0.86580 -0.78875 -0.05764  0.31464 -0.19754 -0.07204
## 
## 
## Site scores (weighted sums of species scores)
## 
##          PC1      PC2       PC3      PC4      PC5      PC6
## 1  -1.283073  1.46362  0.329969  1.17010  0.42054  0.64811
## 2  -1.034856  1.29806  0.047297 -1.40897  0.23982  1.05835
## 3  -0.882710  1.69037 -0.205471  0.53435 -0.97579 -0.24031
## 4  -0.994400  0.94315  0.505607  0.14215 -0.25215 -0.60337
## 5  -2.431354  0.65908 -1.898250 -0.84073 -0.96461 -1.22401
## 7   0.625497  0.59533  0.335832  0.26999  0.46695 -1.53700
## 8   1.479591  0.71474 -2.670917 -1.30598 -0.33714  0.79142
## 9  -0.445864  0.46757 -0.662112  2.69914 -0.97712  0.59770
## 10 -0.761641  0.42997  0.078779 -1.25501  1.05851  0.23628
## 11  0.236561  0.72111 -0.598763  0.65145  1.19091  0.35206
## 12  1.141965  0.47600 -0.093490  0.25066 -0.74183 -0.50609
## 13 -0.266793 -0.81652 -0.519182 -0.08502  0.98419  0.15946
## 14  1.856240 -1.27798 -2.134725  0.31427  0.53622 -1.01871
## 15 -0.551353 -1.32214 -0.026774 -0.19529 -0.31291  0.30059
## 16 -1.000529 -0.09745 -0.165024  0.82015  0.76709 -1.46104
## 17 -1.137819 -0.55812  0.102472  0.12211  0.28418 -0.34966
## 18 -0.635730 -0.10995 -0.083433  0.08852  0.08615  1.19206
## 19 -0.350048 -0.61812  0.238228 -0.74534 -0.27547  0.09870
## 20  0.043308 -0.44619 -0.616293  0.21443 -1.00207  1.03184
## 21 -0.506526 -0.74355  0.520934 -0.32837 -0.52971 -0.12075
## 22 -0.191623 -0.73919  0.679868 -0.07152 -0.12217  2.02785
## 23 -0.542921 -0.42120  0.205271 -0.34908  0.08852 -0.25500
## 24 -0.891549 -1.03600  0.087247 -0.06393  0.27935 -0.67049
## 25 -0.046249 -0.83270  0.886679 -1.66658 -1.29329 -1.16661
## 26  0.287167 -1.00938  0.075905 -0.38458  0.26442  0.03373
## 27  0.062323 -0.58078 -0.420579  0.55599 -0.74345  0.20414
## 28  0.114344 -1.02675 -0.128087 -0.11522  0.48227 -0.18664
## 29  0.131909 -0.90355  0.092727  0.14660  0.57868  0.27300
## 30 -0.242136 -0.88933 -0.900465  0.49281 -0.20915  0.35172
## 31  0.598040  0.66619  0.632611  0.63693 -0.11075  1.10934
## 32  0.291300  0.63471  0.038214  0.41094  0.04070 -0.21870
## 33  0.143544 -0.10039  0.902594  1.70768 -0.31291 -1.76226
## 34 -0.123959  0.49543  0.498420 -1.42250  0.43705 -0.16622
## 35 -0.007017 -0.44478  0.598281 -0.30532  0.06941 -0.01775
## 36  0.325450  0.81900 -0.007016 -0.05172 -0.22195 -0.04243
## 37  0.792494  0.37487 -0.251334  0.32512 -0.90002  0.48860
## 38  0.578051  0.34156  0.563571 -0.15707  0.78511 -0.50635
## 39  0.678964  0.74658  0.548865  0.21253  1.15356  0.84082
## 40  0.552941  0.42893  0.518675 -0.06526  1.10748 -0.31764
## 41  1.096463  0.22508  1.878217 -0.98915 -1.08047 -0.53829
## 42  1.122814  0.65436 -0.062654  0.16617 -2.25180 -0.79297
## 43  0.388900  0.62223  0.979109 -0.59038  0.67362  0.28388
## 44  0.936463  0.08358  0.377874  0.30759  0.86156 -0.93175
## 45  0.630656  1.52392 -1.054410 -1.08937  0.97352 -0.21177
## 46 -0.460918 -0.86694 -0.540394 -0.12777  0.13635 -0.42237
## 47 -0.135065 -0.19959 -0.178417 -0.17257  0.24690 -0.55419
## 48  0.363278  0.39923  0.192034  0.62466  0.65023 -0.15452
## 49 -0.146583 -0.07110 -0.609688 -0.31289 -0.43489  0.59717
## 50 -0.184906 -0.65088 -0.098466  0.17505 -0.44824  0.50493
## 51 -0.329674 -0.27009  0.717715  0.25504 -0.12641  0.24305
## 52  0.668426 -0.10346  0.834806 -0.37667 -1.70236  1.04501
## 53  0.230527 -0.56106  0.072573  0.51022  0.87297  0.49728
## 54  0.208081 -0.77747  0.385568  0.67163  0.59043  1.00980

layout(matrix(1:2, nrow=1))
biplot(pca2,scaling=1)
biplot(pca2,scaling=2)

# note different scaling factors, but solution remains same


##############################
# some follow-up suggestions #
# test PCA-axes for effect of landuse or stream size (as log(Q)) using ANOVA or ANCOVA
# correlate PCA-axes with other potential "controlling" variables (e.g. canopy cover) to give "meta-dimensions" more meaning
# useful function envfit() to relate additional variables to the ordination space

Redundancy analysis (RDA) in R

The classical ordination methods always target dimension reduction

Redundancy analysis (RDA) in R

Two involved matrices: one dependent, one independent.

(Note: For correlation of two matrices see CCorA (Canonical correlation analysis) \(\neq\) CCA.)

Redundancy: The proportion of total variance of in the response variables that can be explained linear combinations of predictors.

Two steps:

  1. MLRs relate each response variable to the independent matrix and predict the response.
  2. The matrix of predicted response variables (same size as original: n objects * p variables) is subject to PCA.

The response variables are constrained to be linear combinations of the predictors first! The PCA can only “ordinate” variation of the responses that is relatable to predictors.

Significance tests for the overall model, for the various RDA-axes and for the individual predictors are available (permutation-based).

Site scores:

Redundancy analysis (RDA) in R

zwc<-scale(wc) # must at least be centered even if dimensionally homogeneous!
xmat<-data.frame(logQ=log(mara.raw$Q),logTDN=log(mara.raw$TDN),canopy=mara.raw$canopy)[-6,]

rda<-rda(zwc~logQ+logTDN+canopy,data=xmat) # take care: confusing X and Y argument names

# actual RDA output check
summary(rda)
## 
## Call:
## rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          18.000       1.00
## Constrained     6.481       0.36
## Unconstrained  11.519       0.64
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                         RDA1    RDA2    RDA3    PC1     PC2     PC3     PC4
## Eigenvalue            5.5244 0.76796 0.18831 4.9837 1.63026 1.19739 0.77385
## Proportion Explained  0.3069 0.04266 0.01046 0.2769 0.09057 0.06652 0.04299
## Cumulative Proportion 0.3069 0.34957 0.36004 0.6369 0.72748 0.79400 0.83699
##                           PC5     PC6     PC7     PC8     PC9     PC10     PC11
## Eigenvalue            0.64230 0.60556 0.46542 0.38373 0.25602 0.153021 0.117674
## Proportion Explained  0.03568 0.03364 0.02586 0.02132 0.01422 0.008501 0.006537
## Cumulative Proportion 0.87268 0.90632 0.93217 0.95349 0.96772 0.976216 0.982754
##                           PC12     PC13     PC14     PC15      PC16      PC17
## Eigenvalue            0.101408 0.079286 0.061438 0.046321 0.0156525 0.0063286
## Proportion Explained  0.005634 0.004405 0.003413 0.002573 0.0008696 0.0003516
## Cumulative Proportion 0.988387 0.992792 0.996205 0.998779 0.9996484 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1   RDA2    RDA3
## Eigenvalue            5.5244 0.7680 0.18831
## Proportion Explained  0.8524 0.1185 0.02906
## Cumulative Proportion 0.8524 0.9709 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  5.531195 
## 
## 
## Species scores
## 
##         RDA1      RDA2      RDA3        PC1        PC2        PC3
## cond -0.6632  0.526208 -0.053197 -7.126e-01 -2.642e-01 -2.366e-01
## pH    0.4359 -0.080604 -0.110311 -2.354e-01  6.035e-01 -2.708e-01
## turb -0.8597  0.004151  0.171732  3.490e-01  4.901e-01 -5.836e-01
## Alk  -0.1085  0.409957 -0.026803 -8.856e-01 -5.159e-01 -4.483e-02
## TSS  -0.7648  0.195626  0.067225  2.585e-01  4.908e-01 -7.306e-01
## PO4  -0.1714 -0.057514 -0.327410  2.628e-01  9.564e-01  3.017e-01
## NH4  -0.6307  0.217307 -0.051937 -4.163e-01 -3.196e-01 -5.981e-01
## NO2  -0.8571 -0.178053  0.281840 -1.528e-01  2.746e-01 -3.129e-01
## NO3  -1.1170 -0.494302 -0.147602  1.649e-01  1.328e-01  3.170e-01
## Br   -0.9585  0.116087 -0.072049 -6.810e-01  2.180e-01  6.036e-02
## Cl   -0.8351 -0.016235  0.006982 -6.587e-01  4.095e-01  4.346e-01
## Fl   -0.7748  0.116617 -0.017050 -7.011e-01  4.096e-01  3.659e-02
## SO4  -0.5605  0.008991  0.140803 -9.589e-01  2.172e-01  3.459e-01
## Na   -0.7506  0.071258  0.039663 -8.902e-01  1.087e-01  6.183e-02
## K    -0.4348  0.085161 -0.092075 -1.085e+00  1.745e-01 -4.157e-02
## Ca   -0.3188  0.430074 -0.088844 -1.078e+00 -7.507e-02 -1.775e-01
## Mg   -0.3587  0.282375 -0.044108 -1.095e+00  4.243e-02 -1.784e-01
## TDN  -1.2252 -0.429552 -0.118224  2.170e-16 -1.802e-16  8.924e-17
## 
## 
## Site scores (weighted sums of species scores)
## 
##         RDA1      RDA2     RDA3      PC1      PC2       PC3
## 1  -0.899757  3.084889 -0.55294 -1.73973 -0.49579  1.405403
## 2  -0.601643  4.018791  0.23349 -1.62902  1.01934 -0.493815
## 3  -0.425217  2.834368  0.10401 -1.40811 -1.10996  0.502236
## 4  -0.756630  1.881205 -1.51654 -1.37687 -0.19517  0.362439
## 5  -2.455504  4.254210  2.42612 -2.05980 -0.84349 -2.485533
## 7   0.861601 -0.365693 -1.72079 -0.14991 -0.56146  0.155296
## 8   1.764345  1.159717  2.64686  0.62914 -1.54289 -2.968104
## 9  -0.482278  0.089585  0.15601 -0.35169 -2.59873  1.153105
## 10 -0.595642  2.398130 -0.19763 -0.94864  1.32755 -0.723328
## 11  0.451013  1.306409 -1.18506 -0.42811 -0.88023 -0.088124
## 12  1.394450 -0.925289  1.78413 -0.03078 -0.38908 -0.245790
## 13 -0.567171 -0.078456  0.53359  0.88444  0.28402  0.167429
## 14  1.469933 -2.003658  3.48847  3.00526 -1.06668 -0.266833
## 15 -1.033131 -1.201150  0.27119  0.59853  0.67023 -0.195910
## 16 -1.140389  0.708121  0.03320 -0.37468 -0.34163  0.496573
## 17 -1.426847  0.273804 -0.49117 -0.56728  0.24291 -0.207985
## 18 -0.726797  0.919409 -0.63159 -0.09927  0.15910  0.396983
## 19 -0.562127 -0.441430 -0.39299  0.45991  0.63986 -0.006655
## 20 -0.146589 -0.580216  1.46528  0.62302 -0.45546 -0.030056
## 21 -0.792760 -0.936524 -1.52046  0.06731  0.43021 -0.279145
## 22 -0.450183 -0.920766 -1.69925  0.35532  0.44301  0.394593
## 23 -0.722434 -0.099725 -1.05913  0.11364  0.32640  0.053559
## 24 -1.294677 -0.621767  0.54340 -0.14052  0.78301 -0.387924
## 25 -0.256636 -1.469854 -1.00794  0.30718  1.42878 -0.901286
## 26 -0.006011 -1.524849  0.44526  0.71389  0.74750 -0.370880
## 27 -0.167762 -1.210237  1.92863  0.73162 -0.47637  0.287856
## 28 -0.214533 -1.385767  0.40575  0.68065  0.48451 -0.293713
## 29 -0.153241 -1.359169  0.41669  0.54181  0.45796  0.113385
## 30 -0.608704 -0.767504  2.40998  0.58176 -0.34287 -0.271359
## 31  0.877051  0.127594 -0.99204 -0.11448 -0.08793  1.149609
## 32  0.524430  0.342174  0.11586 -0.19013 -0.41330  0.388517
## 33  0.100447 -1.912077 -0.27986 -0.07708 -0.63553  1.401936
## 34  0.113385  1.339155 -0.76471 -0.15951  1.13570 -0.158967
## 35 -0.144093 -0.913297 -2.68176  0.14873  0.34036 -0.116528
## 36  0.635959  0.706963 -0.20569 -0.45058 -0.27741 -0.210309
## 37  0.964039 -0.502987  1.56066  0.58492 -0.58515  0.246218
## 38  0.780714 -0.223068 -1.48202 -0.02480  0.36112  0.196388
## 39  1.002097  0.320595 -2.44794 -0.11533 -0.06812  0.544981
## 40  0.777442 -0.009334 -1.95549 -0.09691  0.23264  0.233320
## 41  1.374032 -1.702035  0.01784 -0.03704  1.50770  0.434674
## 42  1.421599 -1.106980  2.45558 -0.07522 -0.69935 -0.519837
## 43  0.700872  0.512436 -1.78976 -0.23677  0.87534  0.442824
## 44  1.074898 -1.074960  0.76215  0.24269  0.26216  0.475054
## 45  1.160824  2.438268 -2.38973  0.14196 -1.14818 -0.702500
## 46 -0.808194 -0.336965  0.64062 -0.06377  0.17315 -1.013755
## 47 -0.197930  0.127281  0.76092 -0.11609  0.31741 -0.241202
## 48  0.528880 -0.060544  0.32122  0.08513 -0.24501  0.789760
## 49 -0.190294  0.489138  1.74631 -0.11103  0.04844 -0.769059
## 50 -0.443365 -0.924883 -0.08718  0.49389 -0.16078  0.116308
## 51 -0.445001 -0.638721 -2.14463  0.13713  0.06390  0.808019
## 52  0.728592 -1.401673  2.01083  0.17276  0.63744  0.369228
## 53  0.045616 -1.090325 -1.06195  0.41535 -0.04511  0.260961
## 54 -0.036682 -1.542340  0.57423  0.45711  0.26592  0.601943
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##         RDA1     RDA2     RDA3      PC1      PC2       PC3
## 1  -0.192647  0.93038  0.92799 -1.73973 -0.49579  1.405403
## 2   0.427478  1.34335  1.22481 -1.62902  1.01934 -0.493815
## 3  -0.146209  0.99448  0.45228 -1.40811 -1.10996  0.502236
## 4  -0.310322 -0.13504  0.14698 -1.37687 -0.19517  0.362439
## 5  -1.484641 -0.21700  0.38344 -2.05980 -0.84349 -2.485533
## 7   0.674710 -0.78138 -0.51445 -0.14991 -0.56146  0.155296
## 8   1.388710 -0.10748 -0.22376  0.62914 -1.54289 -2.968104
## 9  -0.858711 -0.86050 -1.00156 -0.35169 -2.59873  1.153105
## 10  0.407941  1.11874 -0.68765 -0.94864  1.32755 -0.723328
## 11  0.703912 -0.23995 -0.63485 -0.42811 -0.88023 -0.088124
## 12  1.393070 -1.09377  1.50118 -0.03078 -0.38908 -0.245790
## 13 -0.864232  1.33162 -0.10131  0.88444  0.28402  0.167429
## 14  0.070930  1.59769  0.43814  3.00526 -1.06668 -0.266833
## 15 -1.113644 -0.02139  0.08487  0.59853  0.67023 -0.195910
## 16 -0.972647  0.47906 -0.44170 -0.37468 -0.34163  0.496573
## 17 -1.040565 -0.76624 -0.16526 -0.56728  0.24291 -0.207985
## 18 -0.629396  0.77943  0.19352 -0.09927  0.15910  0.396983
## 19 -0.880788  0.52142  0.58591  0.45991  0.63986 -0.006655
## 20 -0.557724  0.29157  0.85510  0.62302 -0.45546 -0.030056
## 21 -0.816524 -0.74571 -1.49309  0.06731  0.43021 -0.279145
## 22 -0.679726 -0.61138 -0.94934  0.35532  0.44301  0.394593
## 23 -0.876393  0.22614  0.17538  0.11364  0.32640  0.053559
## 24 -0.835658 -0.46318 -0.11007 -0.14052  0.78301 -0.387924
## 25 -0.306131 -0.56840 -0.71753  0.30718  1.42878 -0.901286
## 26 -0.050644 -0.33280 -0.12269  0.71389  0.74750 -0.370880
## 27 -0.663953  0.25632  0.79512  0.73162 -0.47637  0.287856
## 28 -0.304155 -0.31610  0.19760  0.68065  0.48451 -0.293713
## 29 -0.229915 -0.50523  0.78395  0.54181  0.45796  0.113385
## 30 -0.750822  0.01951  0.79961  0.58176 -0.34287 -0.271359
## 31  0.919784  0.24149 -0.12906 -0.11448 -0.08793  1.149609
## 32  0.500061  0.28682  0.07581 -0.19013 -0.41330  0.388517
## 33 -0.169995 -1.32278  0.05458 -0.07708 -0.63553  1.401936
## 34  0.327650  1.71656 -0.89375 -0.15951  1.13570 -0.158967
## 35 -0.278900 -0.83339 -1.12100  0.14873  0.34036 -0.116528
## 36  0.827677  0.04073 -0.62665 -0.45058 -0.27741 -0.210309
## 37  0.573226  0.88396 -1.00915  0.58492 -0.58515  0.246218
## 38  0.850153 -0.05214 -0.58813 -0.02480  0.36112  0.196388
## 39  1.044336  0.18841 -1.50261 -0.11533 -0.06812  0.544981
## 40  0.840242 -0.12799 -0.57504 -0.09691  0.23264  0.233320
## 41  1.451280 -0.87165  0.50304 -0.03704  1.50770  0.434674
## 42  1.324854 -1.06452 -0.01889 -0.07522 -0.69935 -0.519837
## 43  0.954237  0.72568 -1.21010 -0.23677  0.87534  0.442824
## 44  1.025094 -0.47834  1.59291  0.24269  0.26216  0.475054
## 45  0.465138  1.60579  0.20294  0.14196 -1.14818 -0.702500
## 46 -0.461071 -1.14745 -0.17978 -0.06377  0.17315 -1.013755
## 47  0.002323 -0.19579  1.34466 -0.11609  0.31741 -0.241202
## 48  0.426895  0.53629 -0.03598  0.08513 -0.24501  0.789760
## 49  0.119534  0.16802 -0.14099 -0.11103  0.04844 -0.769059
## 50 -0.832201 -0.18606  0.31023  0.49389 -0.16078  0.116308
## 51 -0.803215 -0.10485 -0.18123  0.13713  0.06390  0.808019
## 52  0.580101 -0.79935  1.87566  0.17276  0.63744  0.369228
## 53 -0.091596 -0.62895 -0.07307  0.41535 -0.04511  0.260961
## 54 -0.096912 -0.70466 -0.05704  0.45711  0.26592  0.601943
## 
## 
## Biplot scores for constraining variables
## 
##           RDA1    RDA2     RDA3 PC1 PC2 PC3
## logQ    0.5174 -0.6315  0.57757   0   0   0
## logTDN -0.9398 -0.3295 -0.09068   0   0   0
## canopy  0.5328 -0.2110 -0.81953   0   0   0
str(rda)
## List of 12
##  $ colsum        : Named num [1:18] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:18] "cond" "pH" "turb" "Alk" ...
##  $ tot.chi       : num 18
##  $ Ybar          : num [1:53, 1:18] 0.248 0.215 0.216 0.178 0.42 ...
##   ..- attr(*, "scaled:center")= Named num [1:18] 8.64e-16 -1.43e-15 -1.05e-18 -7.64e-16 -1.31e-16 ...
##   .. ..- attr(*, "names")= chr [1:18] "cond" "pH" "turb" "Alk" ...
##   ..- attr(*, "scaled:scale")= Named num [1:18] 0.458 0.269 3.894 0.578 0.82 ...
##   .. ..- attr(*, "names")= chr [1:18] "cond" "pH" "turb" "Alk" ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:53] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:18] "cond" "pH" "turb" "Alk" ...
##   ..- attr(*, "METHOD")= chr "PCA"
##  $ method        : chr "rda"
##  $ call          : language rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##  $ pCCA          : NULL
##  $ CCA           :List of 12
##   ..$ eig      : Named num [1:3] 5.524 0.768 0.188
##   .. ..- attr(*, "names")= chr [1:3] "RDA1" "RDA2" "RDA3"
##   ..$ poseig   : NULL
##   ..$ u        : num [1:53, 1:3] -0.0348 0.0773 -0.0264 -0.0561 -0.2684 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:53] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:3] "RDA1" "RDA2" "RDA3"
##   ..$ v        : num [1:18, 1:3] -0.2164 0.1423 -0.2806 -0.0354 -0.2496 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:18] "cond" "pH" "turb" "Alk" ...
##   .. .. ..$ : chr [1:3] "RDA1" "RDA2" "RDA3"
##   ..$ wa       : num [1:53, 1:3] -0.1627 -0.1088 -0.0769 -0.1368 -0.4439 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:53] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:3] "RDA1" "RDA2" "RDA3"
##   ..$ alias    : NULL
##   ..$ biplot   : num [1:3, 1:3] 0.517 -0.94 0.533 -0.631 -0.329 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:3] "logQ" "logTDN" "canopy"
##   .. .. ..$ : chr [1:3] "RDA1" "RDA2" "RDA3"
##   ..$ rank     : int 3
##   ..$ qrank    : int 3
##   ..$ tot.chi  : num 6.48
##   ..$ QR       :List of 4
##   .. ..$ qr   : num [1:53, 1:3] 14.67002 -0.01451 0.07998 -0.00174 0.07405 ...
##   .. .. ..- attr(*, "assign")= int [1:3] 1 2 3
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:53] "1" "2" "3" "4" ...
##   .. .. .. ..$ : chr [1:3] "logQ" "logTDN" "canopy"
##   .. .. ..- attr(*, "scaled:center")= Named num [1:3] 1.674 -9.205 0.464
##   .. .. .. ..- attr(*, "names")= chr [1:3] "logQ" "logTDN" "canopy"
##   .. ..$ rank : int 3
##   .. ..$ qraux: num [1:3] 1.03 1.18 1.13
##   .. ..$ pivot: int [1:3] 1 2 3
##   .. ..- attr(*, "class")= chr "qr"
##   ..$ envcentre: Named num [1:3] 1.674 -9.205 0.464
##   .. ..- attr(*, "names")= chr [1:3] "logQ" "logTDN" "canopy"
##  $ CA            :List of 6
##   ..$ eig    : Named num [1:17] 4.984 1.63 1.197 0.774 0.642 ...
##   .. ..- attr(*, "names")= chr [1:17] "PC1" "PC2" "PC3" "PC4" ...
##   ..$ poseig : NULL
##   ..$ u      : num [1:53, 1:17] -0.315 -0.295 -0.255 -0.249 -0.372 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:53] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:17] "PC1" "PC2" "PC3" "PC4" ...
##   ..$ v      : num [1:18, 1:17] -0.2449 -0.0809 0.1199 -0.3043 0.0888 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:18] "cond" "pH" "turb" "Alk" ...
##   .. .. ..$ : chr [1:17] "PC1" "PC2" "PC3" "PC4" ...
##   ..$ rank   : int 17
##   ..$ tot.chi: num 11.5
##  $ inertia       : chr "variance"
##  $ regularization: chr "this is a vegan::rda result object"
##  $ terms         :Classes 'terms', 'formula'  language zwc ~ logQ + logTDN + canopy
##   .. ..- attr(*, "variables")= language list(zwc, logQ, logTDN, canopy)
##   .. ..- attr(*, "factors")= int [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:4] "zwc" "logQ" "logTDN" "canopy"
##   .. .. .. ..$ : chr [1:3] "logQ" "logTDN" "canopy"
##   .. ..- attr(*, "term.labels")= chr [1:3] "logQ" "logTDN" "canopy"
##   .. ..- attr(*, "order")= int [1:3] 1 1 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  $ terminfo      :List of 4
##   ..$ terms  :Classes 'terms', 'formula'  language ~logQ + logTDN + canopy
##   .. .. ..- attr(*, "variables")= language list(logQ, logTDN, canopy)
##   .. .. ..- attr(*, "factors")= int [1:3, 1:3] 1 0 0 0 1 0 0 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:3] "logQ" "logTDN" "canopy"
##   .. .. .. .. ..$ : chr [1:3] "logQ" "logTDN" "canopy"
##   .. .. ..- attr(*, "term.labels")= chr [1:3] "logQ" "logTDN" "canopy"
##   .. .. ..- attr(*, "order")= int [1:3] 1 1 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= num 0
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   ..$ assign : int [1:3] 1 2 3
##   ..$ xlev   : Named list()
##   ..$ ordered: Named logi [1:3] FALSE FALSE FALSE
##   .. ..- attr(*, "names")= chr [1:3] "logQ" "logTDN" "canopy"
##  - attr(*, "class")= chr [1:2] "rda" "cca"
?cca.object
## starte den http Server für die Hilfe fertig
RsquareAdj(rda) # redundancy statistic (fractional amount of variation of the response data matrix explained by constraints)
## $r.squared
## [1] 0.3600362
## 
## $adj.r.squared
## [1] 0.3208547

# hypothesis tests #
# testing the first axis (global test)
anova(rda)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance      F Pr(>F)    
## Model     3   6.4807 9.1889  0.001 ***
## Residual 49  11.5193                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda,first=TRUE)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance      F Pr(>F)    
## RDA1      1   5.5244 23.499  0.001 ***
## Residual 49  11.5193                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing all axes sequentially (preceding axes are taken as constraints)
anova(rda,by="axis",model="direct",perm.max=9999,step=1000)
## Permutation test for rda under direct model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## RDA1      1   5.5244 23.4991  0.001 ***
## RDA2      1   0.7680  3.2667  0.050 *  
## RDA3      1   0.1883  0.8010  0.518    
## Residual 49  11.5193                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing the individual terms=constraints
anova(rda,by="terms",model="direct",perm.max=9999,step=1000)  # tests terms sequentially, order matters!
## Permutation test for rda under direct model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## logQ      1   1.8477  7.8596  0.001 ***
## logTDN    1   3.9175 16.6639  0.001 ***
## canopy    1   0.7155  3.0434  0.021 *  
## Residual 49  11.5193                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(rda,by="margin",model="direct",perm.max=9999,step=1000) # tests each term in full model (like drop1() function)
## Permutation test for rda under direct model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = zwc ~ logQ + logTDN + canopy, data = xmat)
##          Df Variance       F Pr(>F)    
## logQ      1   1.0182  4.3312  0.007 ** 
## logTDN    1   2.6904 11.4441  0.001 ***
## canopy    1   0.7155  3.0434  0.027 *  
## Residual 49  11.5193                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

###################
# making triplots #

# again various types of scaling for the plotting step:
# scaling 1 "distance triplot"
# only angles between constraints and responses reflect their correlations (not angles among responses)
# distances among sites reflect their Euclidean distances

# scaling 2 "correlation triplot"
# all angles between constraints and responses reflect correlations 
# distances among sites do not reflect their Euclidean distances

# in both scaling types sites can be projected on constraints and on responses
# factor constraints are shown as centroids instead of arrows, projecting works identical 

# scaling 3 is compromise

# build an RDA scaling type 1 triplot
plot(rda,scaling=1)


(sites<-scores(rda,choices=c(1,2),display="sites",scaling=1)) 
##            RDA1         RDA2
## 1  -0.498460582  0.637194072
## 2  -0.333307201  0.830094727
## 3  -0.235567821  0.585448173
## 4  -0.419169104  0.388569136
## 5  -1.360336512  0.878721129
## 7   0.477322538 -0.075535075
## 8   0.977437916  0.239543485
## 9  -0.267179489  0.018503982
## 10 -0.329982722  0.495341683
## 11  0.249858739  0.269843172
## 12  0.772517930 -0.191121436
## 13 -0.314209914 -0.016205304
## 14  0.814335423 -0.413862255
## 15 -0.572349386 -0.248101430
## 16 -0.631769326  0.146264714
## 17 -0.790465648  0.056555154
## 18 -0.402641922  0.189907030
## 19 -0.311415221 -0.091178758
## 20 -0.081209267 -0.119845614
## 21 -0.439185103 -0.193442177
## 22 -0.249398850 -0.190187243
## 23 -0.400224485 -0.020598570
## 24 -0.717244600 -0.128427989
## 25 -0.142174974 -0.303603256
## 26 -0.003329834 -0.314962678
## 27 -0.092939197 -0.249978470
## 28 -0.118850427 -0.286234700
## 29 -0.084894890 -0.280740961
## 30 -0.337218829 -0.158530543
## 31  0.485881970  0.026354971
## 32  0.290531318  0.070677146
## 33  0.055647198 -0.394945899
## 34  0.062814787  0.276606892
## 35 -0.079826826 -0.188644454
## 36  0.352318114  0.146025624
## 37  0.534072469 -0.103893687
## 38  0.432511681 -0.046075409
## 39  0.555156725  0.066219970
## 40  0.430698723 -0.001927959
## 41  0.761206836 -0.351560915
## 42  0.787558621 -0.228650302
## 43  0.388279575  0.105845304
## 44  0.595488024 -0.222036636
## 45  0.643090631  0.503632310
## 46 -0.447735520 -0.069601236
## 47 -0.109651939  0.026290352
## 48  0.292996889 -0.012505478
## 49 -0.105422063  0.101033121
## 50 -0.245621890 -0.191037745
## 51 -0.246528053 -0.131929977
## 52  0.403636317 -0.289520170
## 53  0.025270757 -0.225210312
## 54 -0.020321586 -0.318575509
## attr(,"const")
## [1] 5.531195

(lcs<-scores(rda,choices=c(1,2),display="lc",scaling=1)) # fitted/constrained site scores
##            RDA1         RDA2
## 1  -0.106725620  0.192172696
## 2   0.236820410  0.277474421
## 3  -0.080999194  0.205412919
## 4  -0.171916633 -0.027892958
## 5  -0.822483293 -0.044821584
## 7   0.373785769 -0.161396377
## 8   0.769338143 -0.022200417
## 9  -0.475721291 -0.177739370
## 10  0.225997452  0.231079381
## 11  0.389963824 -0.049561659
## 12  0.771753557 -0.225921755
## 13 -0.478779804  0.275051081
## 14  0.039294802  0.330009134
## 15 -0.616953195 -0.004418402
## 16 -0.538841264  0.098951894
## 17 -0.576467626 -0.158269346
## 18 -0.348682347  0.160994450
## 19 -0.487951895  0.107700881
## 20 -0.308976357  0.060225476
## 21 -0.452349895 -0.154028956
## 22 -0.376564556 -0.126282122
## 23 -0.485517314  0.046710238
## 24 -0.462950388 -0.095671797
## 25 -0.169594740 -0.117405573
## 26 -0.028056760 -0.068740794
## 27 -0.367826399  0.052943956
## 28 -0.168500544 -0.065292121
## 29 -0.127371851 -0.104357404
## 30 -0.415951385  0.004030800
## 31  0.509555593  0.049880617
## 32  0.277030998  0.059243495
## 33 -0.094176545 -0.273225329
## 34  0.181516631  0.354560391
## 35 -0.154509114 -0.172139937
## 36  0.458528832  0.008412040
## 37  0.317564066  0.182584970
## 38  0.470980327 -0.010770442
## 39  0.578556877  0.038917612
## 40  0.465489775 -0.026436288
## 41  0.804001463 -0.180041364
## 42  0.733962381 -0.219879621
## 43  0.528642419  0.149890940
## 44  0.567896692 -0.098802741
## 45  0.257683843  0.331681637
## 46 -0.255430728 -0.237010393
## 47  0.001287162 -0.040440796
## 48  0.236497729  0.110772850
## 49  0.066221216  0.034706118
## 50 -0.461034789 -0.038432073
## 51 -0.444977082 -0.021657745
## 52  0.321373087 -0.165108083
## 53 -0.050743664 -0.129911681
## 54 -0.053688775 -0.145550870
## attr(,"const")
## [1] 5.531195

(species<-scores(rda,choices=c(1,2),display="sp",scaling=1)*0.5)
##             RDA1        RDA2
## cond -0.59858199  1.27378225
## pH    0.39342255 -0.19511683
## turb -0.77595131  0.01004908
## Alk  -0.09788046  0.99237637
## TSS  -0.69028937  0.47354883
## PO4  -0.15467198 -0.13922271
## NH4  -0.56926510  0.52603200
## NO2  -0.77353865 -0.43101017
## NO3  -1.00809786 -1.19654781
## Br   -0.86504284  0.28100954
## Cl   -0.75374148 -0.03929891
## Fl   -0.69926370  0.28229280
## SO4  -0.50586602  0.02176366
## Na   -0.67742612  0.17249327
## K    -0.39244249  0.20614834
## Ca   -0.28777201  1.04107246
## Mg   -0.32375066  0.68353998
## TDN  -1.10581086 -1.03980988
## attr(,"const")
## [1] 5.531195

(constraints<-scores(rda,choices=c(1,2),display="bp",scaling=1)*2)
##              RDA1       RDA2
## logQ    0.5732272 -0.2608649
## logTDN -1.0412848 -0.1361117
## canopy  0.5903155 -0.0871586
## attr(,"const")
## [1] 5.531195

plot(sites,asp=1,pch=21,bg=landuse,ylim=c(-1.5,1.5))
Arrows(x0=0,y0=0,x1=constraints[,1],y1=constraints[,2],lwd=1.5,col="blue")
text(constraints[,1:2]*1.1,label=rownames(constraints),pos=4,cex=0.8,col="blue")

Arrows(x0=0,y0=0,x1=species[,1],y1=species[,2],lwd=1,arr.length=0)
text(species[,1:2]*1.1,label=rownames(species),pos=4,cex=0.6)

Correspondence Analysis (CA)

Based on one matrix (usually a species or community matrix).

Considers unimodal responses to (unknown) environmental variables.

An indirect GA, resulting gradients are synthetic environmental gradients.

The basis for CA is weighted averaging from environmental and species tables. If env exists, then this can be done to extract bioindicatory information:

\[ u^*=\frac{y_1x_1+y_2x_2+...+y_nx_n}{y_1+y_2+...+y_n} \] A species optimum \(u^*\) is computed as an abundance-weighted means of a specific environmental variable over all sites at which a specific species is present.

This approach works best when:

CA uses a two-way weighted averaging with a theoretical environmental variable iteratively in several steps:

  1. Take arbitrary site scores.
  2. Derive species scores by weighted average of sites scores, for species k (of m): \[ u_k=\sum_{i=1}^n{y_{ki}x_i}/\sum_{i=1}^n{y_{ki}} \]
  3. From the species scores new site scores can be derived, for site i (of n): \[ x_i=\sum_{k=1}^m{y_{ki}u_k}/\sum_{k=1}^n{y_{ki}} \]
  4. Rescaling (standardization) of site and species scores.
  5. Repeat 2-3 several times until stabilisation of site and species scores = first CA axis.
  6. Similar procedure to construct second CA axis (uncorrelated to first).

Canonical Correspondence Analysis (CCA)

Two involved matrices: one dependent, one independent.

In the reciprocal averaging of CA a constraint is included:

The result are axes which inform about species-site relationships, but which also have maximized correlation with linear combinations of (environmental) predictors.

Site scores:

data(varespec) # a R dataset on vegetation
data(varechem) # soil chemistry
head(varespec)
##    Callvulg Empenigr Rhodtome Vaccmyrt Vaccviti Pinusylv Descflex Betupube
## 18     0.55    11.13     0.00     0.00    17.80     0.07     0.00        0
## 15     0.67     0.17     0.00     0.35    12.13     0.12     0.00        0
## 24     0.10     1.55     0.00     0.00    13.47     0.25     0.00        0
## 27     0.00    15.13     2.42     5.92    15.97     0.00     3.70        0
## 23     0.00    12.68     0.00     0.00    23.73     0.03     0.00        0
## 19     0.00     8.92     0.00     2.42    10.28     0.12     0.02        0
##    Vacculig Diphcomp Dicrsp Dicrfusc Dicrpoly Hylosple Pleuschr Polypili
## 18     1.60     2.07   0.00     1.62     0.00      0.0     4.67     0.02
## 15     0.00     0.00   0.33    10.92     0.02      0.0    37.75     0.02
## 24     0.00     0.00  23.43     0.00     1.68      0.0    32.92     0.00
## 27     1.12     0.00   0.00     3.63     0.00      6.7    58.07     0.00
## 23     0.00     0.00   0.00     3.42     0.02      0.0    19.42     0.02
## 19     0.00     0.00   0.00     0.32     0.02      0.0    21.03     0.02
##    Polyjuni Polycomm Pohlnuta Ptilcili Barbhatc Cladarbu Cladrang Cladstel
## 18     0.13     0.00     0.13     0.12     0.00    21.73    21.47     3.50
## 15     0.23     0.00     0.03     0.02     0.00    12.05     8.13     0.18
## 24     0.23     0.00     0.32     0.03     0.00     3.58     5.52     0.07
## 27     0.00     0.13     0.02     0.08     0.08     1.42     7.63     2.55
## 23     2.12     0.00     0.17     1.80     0.02     9.08     9.22     0.05
## 19     1.58     0.18     0.07     0.27     0.02     7.23     4.95    22.08
##    Cladunci Cladcocc Cladcorn Cladgrac Cladfimb Cladcris Cladchlo Cladbotr
## 18     0.30     0.18     0.23     0.25     0.25     0.23     0.00     0.00
## 15     2.65     0.13     0.18     0.23     0.25     1.23     0.00     0.00
## 24     8.93     0.00     0.20     0.48     0.00     0.07     0.10     0.02
## 27     0.15     0.00     0.38     0.12     0.10     0.03     0.00     0.02
## 23     0.73     0.08     1.42     0.50     0.17     1.78     0.05     0.05
## 19     0.25     0.10     0.25     0.18     0.10     0.12     0.05     0.02
##    Cladamau Cladsp Cetreric Cetrisla Flavniva Nepharct Stersp Peltapht Icmaeric
## 18     0.08   0.02     0.02     0.00     0.12     0.02   0.62     0.02        0
## 15     0.00   0.00     0.15     0.03     0.00     0.00   0.85     0.00        0
## 24     0.00   0.00     0.78     0.12     0.00     0.00   0.03     0.00        0
## 27     0.00   0.02     0.00     0.00     0.00     0.00   0.00     0.07        0
## 23     0.00   0.00     0.00     0.00     0.02     0.00   1.58     0.33        0
## 19     0.00   0.00     0.00     0.00     0.02     0.00   0.28     0.00        0
##    Cladcerv Claddefo Cladphyl
## 18        0     0.25        0
## 15        0     1.00        0
## 24        0     0.33        0
## 27        0     0.15        0
## 23        0     1.97        0
## 19        0     0.37        0
head(varechem)
##       N    P     K    Ca    Mg    S    Al   Fe    Mn   Zn  Mo Baresoil Humdepth
## 18 19.8 42.1 139.9 519.4  90.0 32.3  39.0 40.9  58.1  4.5 0.3     43.9      2.2
## 15 13.4 39.1 167.3 356.7  70.7 35.2  88.1 39.0  52.4  5.4 0.3     23.6      2.2
## 24 20.2 67.7 207.1 973.3 209.1 58.1 138.0 35.4  32.1 16.8 0.8     21.2      2.0
## 27 20.6 60.8 233.7 834.0 127.2 40.7  15.4  4.4 132.0 10.7 0.2     18.7      2.9
## 23 23.8 54.5 180.6 777.0 125.8 39.5  24.2  3.0  50.1  6.6 0.3     46.0      3.0
## 19 22.8 40.9 171.4 691.8 151.4 40.8 104.8 17.6  43.6  9.1 0.4     40.5      3.8
##     pH
## 18 2.7
## 15 2.8
## 24 3.0
## 27 2.8
## 23 2.7
## 19 2.7
apply(varespec,1,sum) # approximate 100 (total cover), "absolute" abundance data
##     18     15     24     27     23     19     22     16     28     13     14 
##  89.20  89.82  94.21 125.61  90.46  81.27 109.76  88.52 110.70 101.89  81.65 
##     20     25      7      5      6      3      4      2      9     12     10 
##  64.11  94.06 103.38  94.77 110.90 106.67  84.83 119.13 122.60 119.80 122.37 
##     11     21 
## 112.84  99.17

# correspondence analysis #
# run a CA just based on the species data (unconstrained!)
vare.ca<-cca(X=varespec) # function also used for CCA, but here only one matrix X is supplied

summary(vare.ca,scaling=1)
## 
## Call:
## cca(X = varespec) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total           2.083          1
## Unconstrained   2.083          1
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.5249 0.3568 0.2344 0.19546 0.17762 0.12156 0.11549
## Proportion Explained  0.2520 0.1713 0.1125 0.09383 0.08526 0.05835 0.05544
## Cumulative Proportion 0.2520 0.4233 0.5358 0.62962 0.71489 0.77324 0.82868
##                           CA8     CA9    CA10    CA11    CA12    CA13     CA14
## Eigenvalue            0.08894 0.07318 0.05752 0.04434 0.02546 0.01710 0.014896
## Proportion Explained  0.04269 0.03513 0.02761 0.02129 0.01222 0.00821 0.007151
## Cumulative Proportion 0.87137 0.90650 0.93411 0.95539 0.96762 0.97583 0.982978
##                           CA15     CA16     CA17     CA18     CA19      CA20
## Eigenvalue            0.010160 0.007830 0.006032 0.004008 0.002865 0.0019275
## Proportion Explained  0.004877 0.003759 0.002896 0.001924 0.001375 0.0009253
## Cumulative Proportion 0.987855 0.991614 0.994510 0.996434 0.997809 0.9987341
##                            CA21      CA22      CA23
## Eigenvalue            0.0018074 0.0005864 0.0002434
## Proportion Explained  0.0008676 0.0002815 0.0001168
## Cumulative Proportion 0.9996017 0.9998832 1.0000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##                 CA1       CA2      CA3       CA4        CA5       CA6
## Callvulg  0.0303167 -1.597460  0.11455 -2.894569  0.1376073  2.291129
## Empenigr  0.0751030  0.379305  0.39303  0.023675  0.8568729 -0.400964
## Rhodtome  1.1052309  1.499299  3.04284  0.120106  3.2324306 -0.283510
## Vaccmyrt  1.4614812  1.622935  2.72375  0.231688  0.4604556  0.712538
## Vaccviti  0.1468014  0.313436  0.14696  0.243505  0.6868371 -0.147815
## Pinusylv -0.4820096  0.588517 -0.36020 -0.127094  0.4064754  0.386604
## Descflex  1.5348239  1.218806  1.87562 -0.001340 -1.3136979 -0.070731
## Betupube  0.6694503  1.951826  3.84017  1.389423  7.5959115 -0.244478
## Vacculig -0.0830789 -1.629259  1.05063  0.802648 -0.3058811 -1.625341
## Diphcomp -0.5446464 -1.037570  0.52282  0.940275  0.3682126 -1.082929
## Dicrsp    1.8120408  0.360290 -4.92082  3.088562  1.3867372  0.157815
## Dicrfusc  1.2704743 -0.562978 -0.39718 -2.929542  0.3848272 -2.408710
## Dicrpoly  0.7248118  1.409347  0.80341  1.915549  4.5674148  1.295447
## Hylosple  2.0062408  1.743883  2.27549  0.928884 -3.7648428  2.254851
## Pleuschr  1.3102086  0.583036 -0.01004  0.137298 -1.1216144  0.200422
## Polypili -0.3805097 -1.243904  0.54593  1.477188 -0.7276341 -0.387641
## Polyjuni  1.0133795  0.099043 -2.24697  1.510641  0.7729714 -3.062378
## Polycomm  0.8468241  1.321773  1.13585  1.140723  2.6836594 -0.605038
## Pohlnuta -0.0136453  0.589290 -0.35542  0.135481  0.9369707  0.397246
## Ptilcili  0.4223631  1.598584  3.43474  1.400065  6.3209491  0.198935
## Barbhatc  0.5018348  2.119334  4.57303  1.693188  8.1101807  0.645995
## Cladarbu -0.1531729 -1.483884  0.20024  0.193680  0.0734141  0.358926
## Cladrang -0.5502561 -1.084008  0.40552  0.724060 -0.3357992 -0.335924
## Cladstel -1.4373146  1.077753 -0.44397 -0.375926 -0.2421525  0.004212
## Cladunci  0.8151727 -1.006186 -1.82587 -1.389523  1.6046713  3.675908
## Cladcocc -0.2133215 -0.584429 -0.21434 -0.567886 -0.0003788 -0.145303
## Cladcorn  0.2631227 -0.177858 -0.44464  0.272422  0.3992282 -0.306738
## Cladgrac  0.1956947 -0.311167 -0.23894  0.379013  0.4933026  0.037581
## Cladfimb  0.0009213 -0.161418  0.18463 -0.435908  0.4831233 -0.143751
## Cladcris  0.3373031 -0.470369 -0.05093 -0.823855  0.7182250  0.636140
## Cladchlo -0.6200021  1.207278  0.21889  0.426447  1.9506082  0.120722
## Cladbotr  0.5647242  1.047333  2.65330  0.907734  4.4946805  1.201655
## Cladamau -0.6598144 -1.512880  0.83251  1.577699 -0.0407227 -1.419139
## Cladsp   -0.8209003  0.476164 -0.49752 -0.998241 -0.2393208  0.390785
## Cetreric  0.2458192 -0.689228 -1.68427 -0.131681  0.7439412  2.374535
## Cetrisla -0.3465221  1.362693  0.85897  0.396752  2.7526968  0.396591
## Flavniva -1.4391907 -0.833589 -0.12919  0.007071 -1.4841375  2.956977
## Nepharct  1.6813309  0.199484 -4.33509  2.229917  0.9561223 -5.472858
## Stersp   -0.5172793 -2.280900  0.99775  2.377013 -0.8892757 -1.441228
## Peltapht  0.4035858 -0.043265  0.04538  0.711040  0.1824679 -0.841227
## Icmaeric  0.0378754 -2.419595  0.72135  0.361302 -0.3736424 -2.092136
## Cladcerv -0.9232858 -0.005233 -1.22058  0.305290 -0.8142627  0.414135
## Claddefo  0.5190399 -0.496632 -0.15271 -0.695927  0.9042143  0.909191
## Cladphyl -1.2836161  1.155872 -0.79912 -0.741170 -0.1608002  0.490526
## 
## 
## Site scores (weighted averages of species scores)
## 
##          CA1      CA2       CA3      CA4        CA5      CA6
## 18 -0.108122 -0.53705  0.229574  0.24412  0.1405624 -0.14253
## 15  0.697118 -0.14441 -0.031788 -0.21743 -0.2738522 -0.08146
## 24  0.987603  0.15042 -1.348447  0.80472  0.3095168  0.46773
## 27  0.851765  0.49901  0.443559  0.12277 -0.4814871  0.07589
## 23  0.359881 -0.05608  0.145813  0.15087  0.2405263 -0.17770
## 19  0.003545  0.37017  0.027760  0.06168 -0.1158930 -0.03413
## 22  0.860732 -0.11504  0.110869 -1.02169  0.0772348 -0.60530
## 16  0.636936 -0.33250  0.001120 -0.79797  0.0130769 -0.54049
## 28  1.279352  0.81557  0.670053  0.23137 -0.8929976  0.41783
## 13 -0.195009 -0.80564  0.117686 -0.58286 -0.0007212  0.53071
## 14  0.528532 -0.70420 -0.517771 -0.86836  0.5713441  0.91671
## 20  0.382866 -0.18686 -0.004789  0.10156  0.0458125  0.21087
## 25  0.990715  0.11967 -1.110040  0.44929  0.1885902 -0.70694
## 7  -0.264704 -1.06013  0.334900  0.45973 -0.0326631 -0.19945
## 5  -0.428410 -1.20765  0.374344  0.74970 -0.2596294 -0.30467
## 6  -0.330534 -0.77498  0.130760  0.22391  0.0632686  0.09060
## 3  -0.899601  0.12075 -0.075742  0.03842 -0.1489585 -0.12031
## 4  -0.770294 -0.35351 -0.033779 -0.01795 -0.3007839  0.44303
## 2  -0.992193  0.50319 -0.157505 -0.07070 -0.1065172 -0.09928
## 9  -0.937173  0.78688 -0.258119 -0.19377 -0.0343535 -0.01259
## 12 -0.726413  0.49163 -0.157235 -0.08698 -0.0105774 -0.02801
## 10 -1.002083  0.71239 -0.236526 -0.18643 -0.0231666 -0.04928
## 11 -0.322647 -0.03871 -0.001297  0.09029 -0.1481448  0.06934
## 21  0.259527  0.80746  1.124258  0.36083  1.5437866  0.07051
# summary(vare.ca,scaling=2)
# again two different types of scaling are possible for biplots

# scaling 1 (distances among sites matter)
# distances among sites approximate their chi^2 distance
# close sites have similar species abundances
# a site, which is near a specific species, has a high contribution of that species 

# scaling 2 (relationships among species matter)
# distances among species approximate their chi^2 distance
# close species have similar abundances across sites
# a species, which is near a specific site, is more likely to be found at that site

plot(vare.ca,scaling=1)

plot(vare.ca,scaling=2)

plot(vare.ca,scaling=3) # a compromise

# for any scaling take care when interpreting species close to origin:
# these are "everywhere" or have optimum right at the origin (i.e., optimum with regard to both axis shown)

# for more controlled plotting
species.scores<-scores(vare.ca,display="species",scaling=2)
site.scores<-scores(vare.ca,display="sites",scaling=2)

plot(site.scores,col="black",pch=21,xlim=c(-2,2),ylim=c(-2,2))
text(species.scores,col="red",label=names(varespec),cex=0.7)

# post-hoc fitting of an environmental variable
names(varechem)
##  [1] "N"        "P"        "K"        "Ca"       "Mg"       "S"       
##  [7] "Al"       "Fe"       "Mn"       "Zn"       "Mo"       "Baresoil"
## [13] "Humdepth" "pH"
(ef<-envfit(vare.ca,varechem[,12:13],permutations=1999))
## 
## ***VECTORS
## 
##               CA1      CA2     r2 Pr(>r)   
## Baresoil  0.97947 -0.20161 0.2533  0.058 . 
## Humdepth  0.91602  0.40112 0.4524  0.003 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1999
plot(ef)


###########################
# canonical correspondence analysis #
vare.cca<-cca(Y=varespec,X=varechem) # note strange terminology of X and Y in vegan (don´t ask)
vare.cca<-cca(varespec~.,varechem) # hypothesis tests need formula interface (don´t ask)

summary(vare.cca,scaling=1)
## 
## Call:
## cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe +      Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.0832      1.000
## Constrained    1.4415      0.692
## Unconstrained  0.6417      0.308
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.16285 0.14213 0.11795 0.08903 0.07029
## Proportion Explained  0.2107 0.1401 0.07817 0.06823 0.05662 0.04274 0.03374
## Cumulative Proportion 0.2107 0.3507 0.42890 0.49713 0.55375 0.59649 0.63023
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.02801 0.01495 0.006382 0.004015 0.003139 0.002955
## Cumulative Proportion 0.65825 0.67319 0.679576 0.683592 0.686730 0.689685
##                          CCA14     CA1     CA2     CA3     CA4     CA5     CA6
## Eigenvalue            0.004733 0.19776 0.14193 0.10117 0.07079 0.05330 0.03330
## Proportion Explained  0.002272 0.09493 0.06813 0.04857 0.03398 0.02559 0.01598
## Cumulative Proportion 0.691958 0.78689 0.85502 0.90359 0.93757 0.96315 0.97914
##                            CA7      CA8      CA9
## Eigenvalue            0.018868 0.015104 0.009488
## Proportion Explained  0.009057 0.007251 0.004554
## Cumulative Proportion 0.988195 0.995446 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3   CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.1628 0.1421 0.11795 0.08903 0.07029
## Proportion Explained  0.3045 0.2024 0.1130 0.0986 0.08183 0.06176 0.04877
## Cumulative Proportion 0.3045 0.5069 0.6198 0.7184 0.80027 0.86203 0.91080
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.04049 0.02160 0.009223 0.005803 0.004536 0.004271
## Cumulative Proportion 0.95128 0.97288 0.982107 0.987910 0.992446 0.996716
##                          CCA14
## Eigenvalue            0.004733
## Proportion Explained  0.003284
## Cumulative Proportion 1.000000
## 
## Scaling 1 for species and site scores
## * Sites are scaled proportional to eigenvalues
## * Species are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##              CCA1     CCA2     CCA3      CCA4     CCA5     CCA6
## Callvulg  0.11374 -1.73246  4.15753  1.844837  3.13741 -1.15626
## Empenigr -0.27373  0.14089  0.09036 -1.134550 -0.40226  0.03525
## Rhodtome -1.59033 -0.11156  0.19187 -2.490433 -0.62292 -1.73616
## Vaccmyrt -1.92827  0.56943  0.75259 -0.244263 -1.65624 -2.05452
## Vaccviti -0.23029  0.22315 -0.13141 -0.960947  0.24441  0.02995
## Pinusylv  0.36674  0.48934  0.55326 -0.726273  0.85051 -0.21227
## Descflex -2.17952  0.50020 -0.40165  1.608949 -1.38617  1.28224
## Betupube -1.07326 -0.41989 -0.20569 -6.388345 -0.62954 -5.60316
## Vacculig  0.77560 -2.19992 -0.93608  0.469588 -2.78966  1.04277
## Diphcomp  0.14991 -1.65300 -1.03898 -1.412059 -0.78833  2.08551
## Dicrsp   -1.28302  0.42862 -4.34136  0.691801  4.43282  1.30777
## Dicrfusc -0.75393 -0.76901  2.04376 -0.684761  0.32655  2.14058
## Dicrpoly -0.79564  0.14903 -2.01239 -3.186680  2.23820 -3.43647
## Hylosple -2.75940  1.46965  0.12345  3.602353 -2.66866 -0.74851
## Pleuschr -1.39625  0.62360 -0.02266  0.817213 -0.19077  0.06281
## Polypili  0.21763 -0.84392 -1.27708 -0.747467 -0.15333  0.16978
## Polyjuni -0.91607  0.38917 -0.87255 -0.891253 -1.78446  1.17847
## Polycomm -1.34974  0.59359 -0.58214 -2.854381 -1.19037 -2.60320
## Pohlnuta -0.01435  0.46778 -0.34834 -0.931562  1.23465 -0.32446
## Ptilcili -0.86964 -0.22650 -0.14520 -5.594842 -0.48392 -5.05263
## Barbhatc -1.04773 -0.42524 -0.29330 -6.830157 -0.50320 -6.88497
## Cladarbu  0.31928 -1.31813 -0.06534  0.138503 -0.11811 -0.26229
## Cladrang  0.57516 -1.14185 -0.60438  0.280955 -0.47617  0.10938
## Cladstel  1.36834  1.29986  0.20555  0.179762 -0.04827  0.09185
## Cladunci -0.34820  0.11797 -0.03422 -1.037583  2.65119 -0.48962
## Cladcocc  0.33121 -0.25213  0.31806 -0.205437  0.09828  0.41903
## Cladcorn -0.34025  0.12974 -0.22432 -0.686052 -0.31883  0.57211
## Cladgrac -0.16429 -0.34432 -0.39566 -0.533216  0.70217 -0.07237
## Cladfimb  0.03022 -0.16993  0.47734 -0.696051 -0.10470 -0.11656
## Cladcris -0.20689  0.02979  1.04812 -1.124294  0.40186 -0.43505
## Cladchlo  0.66964  1.02386 -0.68975 -1.528620  0.49217 -0.75368
## Cladbotr -1.02718 -0.35199  0.48348 -3.528217  0.63524 -4.23041
## Cladamau -0.02415 -2.15363 -1.80591 -1.323302 -1.02050  2.39498
## Cladsp    1.03577  0.72454  0.76099  0.741439  1.75911  0.41843
## Cetreric  0.09754 -0.07199 -1.05941  0.315233  2.75328 -0.58261
## Cetrisla  0.24027  0.64936 -0.12182 -2.346145  0.48511 -2.31098
## Flavniva  1.31684 -1.19677 -1.15320  5.202080  1.07346 -7.81575
## Nepharct -1.15139  0.36798 -1.38414 -0.153781 -3.31081  2.49381
## Stersp    0.18370 -2.37391 -2.38790 -0.009846 -1.07525  1.39790
## Peltapht -0.60047  0.31181  0.12300 -0.899163 -0.76856  0.65021
## Icmaeric  0.26085 -2.83828 -1.06550 -0.409685 -1.20472  1.06913
## Cladcerv  1.06877 -0.10889 -0.78377  3.250753  0.01418 -3.50019
## Claddefo -0.45498 -0.03870  0.60324 -1.497542  0.85219 -0.63272
## Cladphyl  1.51291  2.08492  0.04117 -0.268421  0.27480  0.48796
## 
## 
## Site scores (weighted averages of species scores)
## 
##        CCA1     CCA2      CCA3     CCA4       CCA5     CCA6
## 18  0.11823 -0.57251 -0.164982 -0.22892 -0.1940177  0.07213
## 15 -0.64276 -0.10649  0.169910  0.11432  0.0521025  0.23988
## 24 -0.84786  0.25736 -1.189184  0.14813  1.3580797  0.22853
## 27 -0.99432  0.35227  0.034639  0.28730 -0.4232940 -0.02911
## 23 -0.39622 -0.09941 -0.054725 -0.43892 -0.1038877  0.02098
## 19 -0.07306  0.38585  0.006695 -0.02930 -0.1896126 -0.02464
## 22 -0.72351 -0.26482  0.855780 -0.16216  0.0893297  0.55882
## 16 -0.50071 -0.42518  0.666714 -0.05991  0.1632133  0.51821
## 28 -1.48534  0.62159  0.100450  0.70953 -0.6209907 -0.35786
## 13  0.26728 -0.79352  0.904036  0.45978  0.6372512 -0.27314
## 14 -0.30232 -0.37451  0.439688 -0.39404  0.9278442  0.04663
## 20 -0.36984 -0.13664 -0.135727 -0.13735  0.0942856  0.03259
## 25 -0.85602  0.13551 -0.587776 -0.01017  0.3304823  0.65496
## 7   0.36936 -1.08951 -0.372699  0.05638 -0.4616061  0.05740
## 5   0.44060 -1.21454 -0.658393  0.16630 -0.4226872  0.15976
## 6   0.39218 -0.69770 -0.189710 -0.03141 -0.0990151 -0.05450
## 3   0.88634  0.21282 -0.085773  0.09809 -0.2111364  0.08974
## 4   0.77344 -0.30247 -0.083929  0.80863  0.1228687 -0.94716
## 2   0.93351  0.60859  0.004559  0.01574 -0.1379707  0.08149
## 9   0.86982  0.91296  0.096369 -0.05063  0.0005494  0.01469
## 12  0.67010  0.58561  0.034417 -0.09231 -0.0424654  0.05488
## 10  0.93439  0.83587  0.093851 -0.06296 -0.0540425  0.05003
## 11  0.30814  0.02923 -0.059108  0.09765 -0.0281526 -0.01160
## 21 -0.47641  0.23201  0.003915 -1.44448 -0.2880137 -1.21174
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        CCA1     CCA2      CCA3     CCA4      CCA5     CCA6
## 18 -0.28028 -0.71553 -0.198602 -0.35622 -0.016645  0.28042
## 15 -0.12605  0.26839  0.183427 -0.19962 -0.026309 -0.23569
## 24 -0.57190  0.13619 -1.113919  0.21486  1.130853  0.07844
## 27 -1.12491  0.26289 -0.227400  0.40474 -0.210924  0.14882
## 23 -0.52704  0.05792  0.103917 -0.34088 -0.098759  0.13090
## 19 -0.44851  0.54086  0.013494 -0.37832 -0.048521 -0.27997
## 22 -0.54244 -0.36270  0.612068 -0.02209  0.194629  0.66118
## 16 -0.09855 -0.62779  0.413117 -0.16871 -0.053130 -0.07504
## 28 -1.37258  0.59298  0.200796  0.71143 -0.478758 -0.19023
## 13  0.10953 -0.73196  1.049988  0.47179  0.604495 -0.16296
## 14 -0.09320  0.10867  0.313802 -0.33147  0.232443 -0.11453
## 20 -0.45423  0.04379 -0.082409 -0.42046  0.381970 -0.22780
## 25 -0.59995  0.15944 -0.222688 -0.02782 -0.388701  0.24252
## 7   0.91721 -1.04185 -0.323016  0.13738 -0.567910 -0.03541
## 5   0.06432 -1.09164 -0.636766  0.01507 -0.151542  0.29546
## 6   0.27735 -0.30740 -0.130894  0.02489 -0.019959  0.10057
## 3   0.63365  0.06729 -0.206032  0.05714 -0.365798 -0.04820
## 4   0.56735 -0.42871 -0.189592  0.87651  0.160886 -0.84789
## 2   1.01789  0.50232  0.038999  0.09780 -0.003433  0.21274
## 9   1.01611  0.86649 -0.006136 -0.04395  0.239962  0.19821
## 12  0.29646  0.12958  0.378874 -0.10628  0.044242  0.11423
## 10  0.73605  0.86077 -0.016804  0.04149 -0.158371  0.07950
## 11  0.39119  0.19766 -0.018368 -0.05333 -0.024357 -0.11579
## 21 -0.45499 -0.12585 -0.070007 -1.04926 -0.070611 -0.65096
## 
## 
## Biplot scores for constraining variables
## 
##              CCA1     CCA2      CCA3      CCA4      CCA5       CCA6
## N        -0.14766 -0.28570  0.002715  0.066860 -0.086965  0.0304388
## P        -0.21110  0.31268 -0.065374  0.180759  0.063227 -0.0363527
## K        -0.24254  0.16634  0.145204  0.180743  0.111771 -0.0586722
## Ca       -0.29655  0.22782 -0.015240  0.037048  0.105769  0.0129929
## Mg       -0.28817  0.18393 -0.057371  0.040679  0.170979 -0.0017181
## S        -0.01594  0.22454  0.059879  0.167562  0.205056 -0.0496189
## Al        0.50996 -0.02564  0.015177  0.147400  0.055261 -0.1004202
## Fe        0.43000 -0.04760 -0.016976  0.099139 -0.023974 -0.0332227
## Mn       -0.47852  0.12132  0.045621  0.109904 -0.047628  0.0538484
## Zn       -0.23723  0.18092 -0.112151  0.130337  0.212656 -0.0003564
## Mo        0.13523 -0.05582 -0.063359  0.122240  0.177367 -0.0935487
## Baresoil -0.35558 -0.13762  0.055249 -0.196250  0.057225 -0.1051508
## Humdepth -0.46157  0.10891  0.109612 -0.051173 -0.001117 -0.0153218
## pH        0.32935  0.04056 -0.131693  0.007887 -0.049995 -0.0176316
summary(vare.cca,scaling=2)
## 
## Call:
## cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe +      Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.0832      1.000
## Constrained    1.4415      0.692
## Unconstrained  0.6417      0.308
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.16285 0.14213 0.11795 0.08903 0.07029
## Proportion Explained  0.2107 0.1401 0.07817 0.06823 0.05662 0.04274 0.03374
## Cumulative Proportion 0.2107 0.3507 0.42890 0.49713 0.55375 0.59649 0.63023
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.02801 0.01495 0.006382 0.004015 0.003139 0.002955
## Cumulative Proportion 0.65825 0.67319 0.679576 0.683592 0.686730 0.689685
##                          CCA14     CA1     CA2     CA3     CA4     CA5     CA6
## Eigenvalue            0.004733 0.19776 0.14193 0.10117 0.07079 0.05330 0.03330
## Proportion Explained  0.002272 0.09493 0.06813 0.04857 0.03398 0.02559 0.01598
## Cumulative Proportion 0.691958 0.78689 0.85502 0.90359 0.93757 0.96315 0.97914
##                            CA7      CA8      CA9
## Eigenvalue            0.018868 0.015104 0.009488
## Proportion Explained  0.009057 0.007251 0.004554
## Cumulative Proportion 0.988195 0.995446 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3   CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4389 0.2918 0.1628 0.1421 0.11795 0.08903 0.07029
## Proportion Explained  0.3045 0.2024 0.1130 0.0986 0.08183 0.06176 0.04877
## Cumulative Proportion 0.3045 0.5069 0.6198 0.7184 0.80027 0.86203 0.91080
##                          CCA8    CCA9    CCA10    CCA11    CCA12    CCA13
## Eigenvalue            0.05836 0.03114 0.013294 0.008364 0.006538 0.006156
## Proportion Explained  0.04049 0.02160 0.009223 0.005803 0.004536 0.004271
## Cumulative Proportion 0.95128 0.97288 0.982107 0.987910 0.992446 0.996716
##                          CCA14
## Eigenvalue            0.004733
## Proportion Explained  0.003284
## Cumulative Proportion 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##               CCA1     CCA2      CCA3      CCA4      CCA5      CCA6
## Callvulg  0.075347 -0.93581  1.677742  0.695507  1.077518 -0.345001
## Empenigr -0.181340  0.07610  0.036462 -0.427727 -0.138153  0.010517
## Rhodtome -1.053549 -0.06026  0.077428 -0.938897 -0.213938 -0.518031
## Vaccmyrt -1.277428  0.30759  0.303704 -0.092088 -0.568820 -0.613023
## Vaccviti -0.152563  0.12054 -0.053031 -0.362279  0.083942  0.008938
## Pinusylv  0.242956  0.26432  0.223265 -0.273806  0.292102 -0.063335
## Descflex -1.443872  0.27019 -0.162082  0.606576 -0.476067  0.382590
## Betupube -0.711004 -0.22681 -0.083007 -2.408417 -0.216212 -1.671857
## Vacculig  0.513817 -1.18831 -0.377748  0.177035 -0.958084  0.311138
## Diphcomp  0.099310 -0.89289 -0.419273 -0.532348 -0.270745  0.622270
## Dicrsp   -0.849964  0.23153 -1.751924  0.260810  1.522412  0.390210
## Dicrfusc -0.499460 -0.41539  0.824743 -0.258156  0.112149  0.638702
## Dicrpoly -0.527090  0.08050 -0.812083 -1.201383  0.768689 -1.025365
## Hylosple -1.828026  0.79385  0.049816  1.358093 -0.916528 -0.223338
## Pleuschr -0.924978  0.33684 -0.009146  0.308091 -0.065518  0.018741
## Polypili  0.144172 -0.45586 -0.515356 -0.281796 -0.052660  0.050659
## Polyjuni -0.606869  0.21021 -0.352109 -0.336004 -0.612858  0.351629
## Polycomm -0.894165  0.32063 -0.234919 -1.076106 -0.408823 -0.776736
## Pohlnuta -0.009508  0.25268 -0.140571 -0.351201  0.424031 -0.096811
## Ptilcili -0.576115 -0.12234 -0.058593 -2.109265 -0.166198 -1.507591
## Barbhatc -0.694092 -0.22970 -0.118360 -2.574980 -0.172821 -2.054320
## Cladarbu  0.211517 -0.71201 -0.026366  0.052216 -0.040564 -0.078262
## Cladrang  0.381030 -0.61678 -0.243893  0.105921 -0.163536  0.032637
## Cladstel  0.906486  0.70213  0.082949  0.067771 -0.016579  0.027407
## Cladunci -0.230671  0.06372 -0.013810 -0.391170  0.910527 -0.146092
## Cladcocc  0.219419 -0.13619  0.128350 -0.077450  0.033754  0.125028
## Cladcorn -0.225404  0.07008 -0.090524 -0.258643 -0.109501  0.170706
## Cladgrac -0.108836 -0.18599 -0.159664 -0.201023  0.241156 -0.021594
## Cladfimb  0.020022 -0.09179  0.192626 -0.262413 -0.035959 -0.034780
## Cladcris -0.137056  0.01609  0.422960 -0.423861  0.138016 -0.129810
## Cladchlo  0.443621  0.55305 -0.278345 -0.576292  0.169030 -0.224882
## Cladbotr -0.680481 -0.19013  0.195105 -1.330144  0.218169 -1.262258
## Cladamau -0.015996 -1.16331 -0.728763 -0.498887 -0.350481  0.714608
## Cladsp    0.686166  0.39137  0.307091  0.279524  0.604150  0.124850
## Cetreric  0.064619 -0.03889 -0.427516  0.118844  0.945590 -0.173838
## Cetrisla  0.159171  0.35076 -0.049161 -0.884501  0.166607 -0.689545
## Flavniva  0.872373 -0.64645 -0.465365  1.961193  0.368671 -2.332045
## Nepharct -0.762768  0.19877 -0.558560 -0.057976 -1.137069  0.744096
## Stersp    0.121697 -1.28229 -0.963619 -0.003712 -0.369284  0.417103
## Peltapht -0.397796  0.16843  0.049634 -0.338986 -0.263955  0.194009
## Icmaeric  0.172805 -1.53313 -0.429975 -0.154452 -0.413750  0.319003
## Cladcerv  0.708032 -0.05882 -0.316283  1.225539  0.004871 -1.044377
## Claddefo -0.301412 -0.02090  0.243431 -0.564576  0.292677 -0.188788
## Cladphyl  1.002262  1.12620  0.016613 -0.101195  0.094379  0.145598
## 
## 
## Site scores (weighted averages of species scores)
## 
##       CCA1     CCA2      CCA3     CCA4     CCA5     CCA6
## 18  0.1785 -1.05988 -0.408835 -0.60721 -0.56492  0.24175
## 15 -0.9702 -0.19714  0.421046  0.30324  0.15171  0.80394
## 24 -1.2798  0.47645 -2.946863  0.39292  3.95433  0.76592
## 27 -1.5009  0.65216  0.085837  0.76207 -1.23251 -0.09756
## 23 -0.5981 -0.18404 -0.135611 -1.16425 -0.30249  0.07033
## 19 -0.1103  0.71431  0.016591 -0.07773 -0.55210 -0.08258
## 22 -1.0921 -0.49026  2.120668 -0.43014  0.26010  1.87287
## 16 -0.7558 -0.78712  1.652152 -0.15892  0.47523  1.73677
## 28 -2.2421  1.15075  0.248921  1.88204 -1.80814 -1.19935
## 13  0.4035 -1.46904  2.240249  1.21956  1.85549 -0.91541
## 14 -0.4563 -0.69333  1.089571 -1.04519  2.70161  0.15628
## 20 -0.5583 -0.25296 -0.336340 -0.36433  0.27453  0.10923
## 25 -1.2922  0.25087 -1.456542 -0.02698  0.96227  2.19508
## 7   0.5576 -2.01700 -0.923568  0.14954 -1.34406  0.19237
## 5   0.6651 -2.24847 -1.631533  0.44110 -1.23074  0.53544
## 6   0.5920 -1.29165 -0.470112 -0.08331 -0.28830 -0.18265
## 3   1.3379  0.39399 -0.212551  0.26020 -0.61477  0.30075
## 4   1.1675 -0.55997 -0.207980  2.14490  0.35776 -3.17436
## 2   1.4091  1.12669  0.011297  0.04175 -0.40173  0.27311
## 9   1.3130  1.69016  0.238808 -0.13429  0.00160  0.04923
## 12  1.0115  1.08413  0.085287 -0.24485 -0.12365  0.18392
## 10  1.4105  1.54744  0.232569 -0.16699 -0.15736  0.16768
## 11  0.4651  0.05411 -0.146473  0.25902 -0.08197 -0.03886
## 21 -0.7191  0.42952  0.009702 -3.83149 -0.83861 -4.06109
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        CCA1     CCA2     CCA3     CCA4      CCA5    CCA6
## 18 -0.42308 -1.32466 -0.49215 -0.94489 -0.048464  0.9398
## 15 -0.19026  0.49687  0.45454 -0.52951 -0.076603 -0.7899
## 24 -0.86328  0.25213 -2.76035  0.56993  3.292710  0.2629
## 27 -1.69805  0.48669 -0.56351  1.07358 -0.614147  0.4988
## 23 -0.79557  0.10723  0.25751 -0.90419 -0.287557  0.4387
## 19 -0.67702  1.00130  0.03344 -1.00351 -0.141279 -0.9383
## 22 -0.81881 -0.67147  1.51674 -0.05858  0.566703  2.2159
## 16 -0.14877 -1.16222  1.02373 -0.44751 -0.154699 -0.2515
## 28 -2.07190  1.09778  0.49758  1.88707 -1.394002 -0.6375
## 13  0.16534 -1.35508  2.60193  1.25142  1.760111 -0.5461
## 14 -0.14069  0.20118  0.77762 -0.87922  0.676806 -0.3838
## 20 -0.68566  0.08107 -0.20421 -1.11529  1.112185 -0.7635
## 25 -0.90562  0.29517 -0.55183 -0.07379 -1.131782  0.8128
## 7   1.38453 -1.92877 -0.80045  0.36440 -1.653585 -0.1187
## 5   0.09709 -2.02095 -1.57794  0.03999 -0.441247  0.9902
## 6   0.41866 -0.56908 -0.32436  0.06603 -0.058116  0.3371
## 3   0.95649  0.12458 -0.51056  0.15157 -1.065096 -0.1616
## 4   0.85641 -0.79366 -0.46982  2.32495  0.468453 -2.8417
## 2   1.53650  0.92994  0.09664  0.25941 -0.009995  0.7130
## 9   1.53381  1.60412 -0.01520 -0.11658  0.698700  0.6643
## 12  0.44751  0.23990  0.93887 -0.28191  0.128819  0.3828
## 10  1.11107  1.59354 -0.04164  0.11005 -0.461130  0.2664
## 11  0.59050  0.36592 -0.04552 -0.14145 -0.070919 -0.3881
## 21 -0.68681 -0.23299 -0.17348 -2.78317 -0.205599 -2.1817
## 
## 
## Biplot scores for constraining variables
## 
##              CCA1     CCA2      CCA3     CCA4      CCA5      CCA6
## N        -0.22290 -0.52891  0.006729  0.17735 -0.253216  0.102014
## P        -0.31866  0.57886 -0.162001  0.47947  0.184099 -0.121835
## K        -0.36612  0.30794  0.359824  0.47942  0.325444 -0.196637
## Ca       -0.44764  0.42176 -0.037765  0.09827  0.307969  0.043545
## Mg       -0.43499  0.34051 -0.142169  0.10790  0.497841 -0.005758
## S        -0.02406  0.41570  0.148384  0.44446  0.597063 -0.166296
## Al        0.76978 -0.04747  0.037610  0.39098  0.160905 -0.336554
## Fe        0.64909 -0.08811 -0.042067  0.26297 -0.069806 -0.111345
## Mn       -0.72232  0.22460  0.113052  0.29152 -0.138680  0.180471
## Zn       -0.35810  0.33493 -0.277916  0.34572  0.619191 -0.001195
## Mo        0.20413 -0.10334 -0.157007  0.32424  0.516439 -0.313525
## Baresoil -0.53675 -0.25477  0.136910 -0.52055  0.166621 -0.352409
## Humdepth -0.69673  0.20163  0.271625 -0.13574 -0.003252 -0.051350
## pH        0.49716  0.07509 -0.326341  0.02092 -0.145569 -0.059091
# again two different types of scaling are possible for triplots

# hypothesis tests #
# testing the first axis (global test)
anova(vare.cca)
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## Model    14   1.44148 1.4441  0.032 *
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vare.cca,first=TRUE)
## Permutation test for cca under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## CCA1      1   0.43887 6.1551  0.041 *
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing all axes sequentially (preceding axes are taken as constraints)
anova(vare.cca,by="axis",model="direct",perm.max=9999,step=1000)
## Permutation test for cca under direct model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## CCA1      1   0.43887 6.1551  0.029 *
## CCA2      1   0.29178 4.0921  0.401  
## CCA3      1   0.16285 2.2839  0.962  
## CCA4      1   0.14213 1.9934  0.970  
## CCA5      1   0.11795 1.6543  0.994  
## CCA6      1   0.08903 1.2486  1.000  
## CCA7      1   0.07029 0.9859  1.000  
## CCA8      1   0.05836 0.8185  1.000  
## CCA9      1   0.03114 0.4367  1.000  
## CCA10     1   0.01329 0.1865  1.000  
## CCA11     1   0.00836 0.1173  1.000  
## CCA12     1   0.00654 0.0917  1.000  
## CCA13     1   0.00616 0.0863  1.000  
## CCA14     1   0.00473 0.0664  1.000  
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# testing the individual terms=constraints
anova(vare.cca,by="terms",model="direct",perm.max=9999,step=1000)  # tests terms sequentially, order matters!
## Permutation test for cca under direct model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## N         1   0.13001 1.8234  0.053 .
## P         1   0.17996 2.5240  0.018 *
## K         1   0.13643 1.9134  0.065 .
## Ca        1   0.08614 1.2081  0.279  
## Mg        1   0.06760 0.9481  0.460  
## S         1   0.17545 2.4606  0.012 *
## Al        1   0.14797 2.0753  0.035 *
## Fe        1   0.04982 0.6987  0.723  
## Mn        1   0.05080 0.7124  0.748  
## Zn        1   0.07167 1.0052  0.452  
## Mo        1   0.09585 1.3442  0.219  
## Baresoil  1   0.08683 1.2178  0.267  
## Humdepth  1   0.09502 1.3327  0.228  
## pH        1   0.06794 0.9529  0.472  
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(vare.cca,by="margin",model="direct",perm.max=9999,step=1000) # tests each term in full model (like drop1() function)
## Permutation test for cca under direct model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = varespec ~ N + P + K + Ca + Mg + S + Al + Fe + Mn + Zn + Mo + Baresoil + Humdepth + pH, data = varechem)
##          Df ChiSquare      F Pr(>F)  
## N         1   0.06302 0.8838  0.548  
## P         1   0.10327 1.4484  0.177  
## K         1   0.10385 1.4565  0.140  
## Ca        1   0.06169 0.8652  0.589  
## Mg        1   0.08830 1.2384  0.258  
## S         1   0.10388 1.4569  0.166  
## Al        1   0.06876 0.9644  0.499  
## Fe        1   0.04765 0.6683  0.748  
## Mn        1   0.08110 1.1375  0.304  
## Zn        1   0.06448 0.9043  0.512  
## Mo        1   0.07761 1.0885  0.351  
## Baresoil  1   0.08966 1.2574  0.272  
## Humdepth  1   0.13548 1.9001  0.050 *
## pH        1   0.06794 0.9529  0.490  
## Residual  9   0.64171                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# quite a lot of variables in the constraining matrix, maybe selection would be adequate
# --> function ordistep()

###################
# making triplots #

# again various types of scaling for the plotting step:
# scaling 1 "distance triplot"
# sites can be projected on constraints
# sites close to centroid of factor constraint are more likely to possess the specific state (factor level)
# distances among sites reflect their Chi^2 distances

# scaling 2 "correlation triplot"
# species can be projected on constraints (to give their optimum)
# species close to centroid of factor constraint are more likely to be found in the respective sites
# distances among sites do not reflect their Chi^2 distances

# scaling 3 is compromise

plot(vare.cca,scaling=1,display=c("species","sites"))
plot(vare.cca,scaling=1,display=c("species","sites","bp"))

plot(vare.cca,scaling=2)

plot(vare.cca,scaling=3) # a compromise

# for any scaling take care when interpreting species close to origin:
# these are "everywhere" or have optimum right at the origin (i.e., optimum with regard to both axis shown)

# for more controlled plotting compute scores individually

(species<-scores(vare.cca,display="species",scaling=2))
##                  CCA1        CCA2
## Callvulg  0.075346828 -0.93580864
## Empenigr -0.181339634  0.07610159
## Rhodtome -1.053549303 -0.06026078
## Vaccmyrt -1.277428382  0.30758571
## Vaccviti -0.152563159  0.12053851
## Pinusylv  0.242955730  0.26432438
## Descflex -1.443872268  0.27018775
## Betupube -0.711004196 -0.22680891
## Vacculig  0.513817276 -1.18831302
## Diphcomp  0.099309627 -0.89289046
## Dicrsp   -0.849963650  0.23152562
## Dicrfusc -0.499459808 -0.41538963
## Dicrpoly -0.527089621  0.08050066
## Hylosple -1.828025570  0.79385044
## Pleuschr -0.924977860  0.33684271
## Polypili  0.144171903 -0.45585592
## Polyjuni -0.606869478  0.21021273
## Polycomm -0.894165298  0.32063347
## Pohlnuta -0.009508326  0.25267888
## Ptilcili -0.576115080 -0.12234488
## Barbhatc -0.694091803 -0.22969682
## Cladarbu  0.211516504 -0.71200604
## Cladrang  0.381029817 -0.61678185
## Cladstel  0.906485693  0.70213449
## Cladunci -0.230670950  0.06372124
## Cladcocc  0.219418695 -0.13619040
## Cladcorn -0.225403730  0.07008025
## Cladgrac -0.108836318 -0.18598627
## Cladfimb  0.020022402 -0.09178764
## Cladcris -0.137056217  0.01608948
## Cladchlo  0.443621232  0.55305019
## Cladbotr -0.680480562 -0.19013376
## Cladamau -0.015995645 -1.16331056
## Cladsp    0.686166357  0.39136838
## Cetreric  0.064619228 -0.03888872
## Cetrisla  0.159171117  0.35076227
## Flavniva  0.872372656 -0.64644751
## Nepharct -0.762767923  0.19876924
## Stersp    0.121696721 -1.28229477
## Peltapht -0.397796035  0.16842583
## Icmaeric  0.172805490 -1.53313439
## Cladcerv  0.708032101 -0.05881737
## Claddefo -0.301411810 -0.02090192
## Cladphyl  1.002261821  1.12619548

(lcs<-scores(vare.cca,display="lc",scaling=2)) # fitted site scores
##           CCA1        CCA2
## 18 -0.42308410 -1.32465581
## 15 -0.19026462  0.49687345
## 24 -0.86328447  0.25212611
## 27 -1.69805066  0.48669129
## 23 -0.79556975  0.10723446
## 19 -0.67702233  1.00129602
## 22 -0.81881442 -0.67146524
## 16 -0.14876518 -1.16222054
## 28 -2.07190438  1.09778134
## 13  0.16533766 -1.35507570
## 14 -0.14068713  0.20117513
## 20 -0.68566132  0.08107181
## 25 -0.90561895  0.29516717
## 7   1.38453057 -1.92876582
## 5   0.09708596 -2.02094655
## 6   0.41865544 -0.56908384
## 3   0.95649196  0.12457945
## 4   0.85641341 -0.79366326
## 2   1.53649728  0.92993610
## 9   1.53380971  1.60412304
## 12  0.44750585  0.23989837
## 10  1.11106621  1.59354430
## 11  0.59049933  0.36592318
## 21 -0.68681020 -0.23299449

(sites<-scores(vare.cca,display="sites",scaling=2)) # unfitted site scores
##          CCA1        CCA2
## 18  0.1784733 -1.05988423
## 15 -0.9702382 -0.19713866
## 24 -1.2798478  0.47644978
## 27 -1.5009195  0.65215591
## 23 -0.5980933 -0.18403623
## 19 -0.1102881  0.71431421
## 22 -1.0921288 -0.49026245
## 16 -0.7558244 -0.78712482
## 28 -2.2421137  1.15075172
## 13  0.4034539 -1.46904155
## 14 -0.4563466 -0.69332916
## 20 -0.5582740 -0.25295922
## 25 -1.2921584  0.25086578
## 7   0.5575516 -2.01700439
## 5   0.6650822 -2.24846553
## 6   0.5919915 -1.29165187
## 3   1.3379214  0.39399435
## 4   1.1674982 -0.55996655
## 2   1.4091299  1.12668910
## 9   1.3129824  1.69015916
## 12  1.0115068  1.08413112
## 10  1.4104517  1.54743675
## 11  0.4651334  0.05410696
## 21 -0.7191331  0.42951771

(constraints<-scores(vare.cca,choices=c(1,2),display="bp",scaling=2))
##                 CCA1        CCA2
## N        -0.22289549 -0.52891359
## P        -0.31865680  0.57886019
## K        -0.36611630  0.30794063
## Ca       -0.44763820  0.42176452
## Mg       -0.43499067  0.34051494
## S        -0.02405594  0.41569739
## Al        0.76977830 -0.04747443
## Fe        0.64908911 -0.08811470
## Mn       -0.72232495  0.22459941
## Zn       -0.35809796  0.33492792
## Mo        0.20413234 -0.10334347
## Baresoil -0.53674724 -0.25477013
## Humdepth -0.69673051  0.20163085
## pH        0.49715734  0.07508792

plot(sites,col="black",pch=21,xlim=c(-2,2),ylim=c(-2,2))
text(species,col="red",label=names(varespec),cex=0.7)
Arrows(x0=0,y0=0,x1=constraints[,1],y1=constraints[,2],lwd=1.5,col="blue")
text(constraints[,1:2]*1.1,label=rownames(constraints),pos=4,cex=0.8,col="blue")